diff --git a/CMakeLists.txt b/CMakeLists.txt index f9b2dd0..7a13543 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2020 Robert Grupp +# Copyright (c) 2020-2021 Robert Grupp # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -22,13 +22,13 @@ cmake_minimum_required(VERSION 3.13.0 FATAL_ERROR) -if (APPLE) - # Target MacOS 10.14, Mojave - set (CMAKE_OSX_DEPLOYMENT_TARGET "10.14" CACHE STRING "MacOS Deployment Target") -endif () +#if (APPLE) +# # Target MacOS 10.14, Mojave +# set (CMAKE_OSX_DEPLOYMENT_TARGET "10.14" CACHE STRING "MacOS Deployment Target") +#endif () # Version format is YYYY.MM.DD for a release and YYYY.MM.DD.1 for devel after the release -project(xreg VERSION 2020.12.13.0) +project(xreg VERSION 2020.12.13.1) set (CMAKE_CXX_STANDARD "11" CACHE STRING "C++ Standard (needs at least 11)") mark_as_advanced(CMAKE_CXX_STANDARD) diff --git a/README.md b/README.md index 6ef23be..feaa211 100644 --- a/README.md +++ b/README.md @@ -70,7 +70,9 @@ The [docker](docker) directory demonstrates how Docker may be used to build the * Highly recomended for GPU acceleration: OpenCL (1.x) * Only needed at runtime on Windows and Linux and is typically provided with your graphics drivers or CUDA SDK * Included with MacOS - * Optional: [ffmpeg](https://ffmpeg.org) is used for writing videos when it is found in the system path. The OpenCV video writer is used if ffmpeg is not available. + * Optional: [ffmpeg](https://ffmpeg.org) is used for writing videos when it is found in the system path. When ffmpeg is not found, the following fallbacks are used: + * MacOS: a video writer using AVFoundation, + * Windows, Linux, something else: a writer using OpenCV. ## Testing Functional testing is available in the form of a [python script](tests/wiki_cmds.py) that runs the commands found on the [wiki walkthrough](https://github.com/rg2/xreg/wiki#walkthrough). diff --git a/apps/hip_surgery/pelvis_single_view_regi_2d_3d/xreg_hip_surg_pelvis_single_view_regi_2d_3d_main.cpp b/apps/hip_surgery/pelvis_single_view_regi_2d_3d/xreg_hip_surg_pelvis_single_view_regi_2d_3d_main.cpp index 0fd871b..92512e4 100644 --- a/apps/hip_surgery/pelvis_single_view_regi_2d_3d/xreg_hip_surg_pelvis_single_view_regi_2d_3d_main.cpp +++ b/apps/hip_surgery/pelvis_single_view_regi_2d_3d/xreg_hip_surg_pelvis_single_view_regi_2d_3d_main.cpp @@ -87,6 +87,10 @@ int main(int argc, char* argv[]) "Do NOT convert RAS to LPS (or LPS to RAS) for the 3D landmarks; " "RAS to LPS conversion negates the first and second components.") << false; + + po.add("no-log-remap", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-log-remap", + "Do NOT perform log remapping of the projection intensities during pre-processing.") + << false; po.add_backend_flags(); @@ -130,6 +134,8 @@ int main(int argc, char* argv[]) const unsigned char pelvis_label = static_cast(po.get("pelvis-label").as_uint32()); + const bool no_log_remap = po.get("no-log-remap"); + ////////////////////////////////////////////////////////////////////////////// // Read in input intensity volume @@ -167,6 +173,7 @@ int main(int argc, char* argv[]) PrintLandmarkMap(lands_3d, vout); ProjPreProc proj_preproc; + proj_preproc.params.no_log_remap = no_log_remap; { vout << "reading projection data..." << std::endl; diff --git a/apps/image_io/CMakeLists.txt b/apps/image_io/CMakeLists.txt index 82131de..557abc5 100644 --- a/apps/image_io/CMakeLists.txt +++ b/apps/image_io/CMakeLists.txt @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2020 Robert Grupp +# Copyright (c) 2020-2021 Robert Grupp # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -28,4 +28,9 @@ add_subdirectory(extract_nii_from_proj_data) add_subdirectory(add_lands_to_proj_data) add_subdirectory(draw_xray_scene) add_subdirectory(regi2d3d_replay) +add_subdirectory(convert_rad_raw_to_proj_data) +add_subdirectory(convert_sta_raw_to_itk) +add_subdirectory(convert_radiograph_dicom_to_proj_data) +add_subdirectory(make_video_from_image_dir) +add_subdirectory(remap_dicom_dir) diff --git a/apps/image_io/add_lands_to_proj_data/xreg_add_lands_to_proj_data_main.cpp b/apps/image_io/add_lands_to_proj_data/xreg_add_lands_to_proj_data_main.cpp index cb88c70..1fea3b1 100644 --- a/apps/image_io/add_lands_to_proj_data/xreg_add_lands_to_proj_data_main.cpp +++ b/apps/image_io/add_lands_to_proj_data/xreg_add_lands_to_proj_data_main.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -55,6 +55,18 @@ int main(int argc, char* argv[]) po.add("no-pat-rot-up", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-pat-rot-up", "Ignore any flags for rotating the image to achieve patient \"up\" orientation.") << false; + + po.add("fcsv-spacing", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_DOUBLE, "fcsv-spacing", + "Default (isotopic) pixel spacing to assume when parsing FCSV landmarks when no 2D pixel " + "spacing is provided by the 2D image metadata.") + << 1.0; + + po.add("use-pd-spacing", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "use-pd-spacing", + "Always use the pixel spacings encoded in the projection data file when interpreting " + "FCSV files, even when the pixel spacings were not explicitly provided by the source " + "image format. This is useful when an FCSV file was created using a NIFTI file which " + "was extracted from a projection data file with estimated spacings.") + << false; try { @@ -78,6 +90,10 @@ int main(int argc, char* argv[]) const bool ignore_pat_rot_up = po.get("no-pat-rot-up"); + const double default_fcsv_spacing = po.get("fcsv-spacing"); + + const bool use_pd_spacing = po.get("use-pd-spacing"); + vout << "opening proj data for reading/writing..." << std::endl; H5::H5File h5(po.pos_args()[0], H5F_ACC_RDWR); @@ -123,6 +139,27 @@ int main(int argc, char* argv[]) CoordScalar spacing_for_phys_to_ind_x = cur_cam.det_col_spacing; CoordScalar spacing_for_phys_to_ind_y = cur_cam.det_row_spacing; + // The following conditional addresses the UNCOMMON problem listed below: + // Suppose that our source data does not have explicit metadata listing pixel spacings, + // so the conversion to projection data HDF5 format used some additional information to + // estimate the pixel spacings. Now, also suppose that 3D Slicer is not able to read in + // the source data - this has happened with certain DICOMs from TCIA with RF modality. + // In order to annotate landmarks, a NIFTI (.nii.gz) file is extracted from the projection + // HDF5 file and then loaded into 3D Slicer for landmark annotation. This NIFTI file has + // the estimated pixel spacings embedded, and therefore the FCSV annotations are physical + // points calculated using these spacings. When the contents of the FCSV are then loaded + // back into the HDF5 projection data (e.g. using THIS TOOL), the current behavior is to use + // a user-provided spacing when the original source data did not have explicit pixel spacings + // provided. Although this is the most common case and default desired behavior, for the + // previously identified scenario we need to use the estimated pixel spacings present in the + // HDF5 file. + if (!use_pd_spacing && cur_proj_meta.det_spacings_from_orig_meta && + !*cur_proj_meta.det_spacings_from_orig_meta) + { + spacing_for_phys_to_ind_x = default_fcsv_spacing; + spacing_for_phys_to_ind_y = default_fcsv_spacing; + } + bool need_to_rot = false; const bool check_rot_field = !ignore_pat_rot_up && cur_proj_meta.rot_to_pat_up; diff --git a/apps/image_io/convert_rad_raw_to_proj_data/CMakeLists.txt b/apps/image_io/convert_rad_raw_to_proj_data/CMakeLists.txt new file mode 100644 index 0000000..1bb9a38 --- /dev/null +++ b/apps/image_io/convert_rad_raw_to_proj_data/CMakeLists.txt @@ -0,0 +1,30 @@ +# MIT License +# +# Copyright (c) 2021 Robert Grupp +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +set(EXE_NAME "${XREG_EXE_PREFIX}convert-rad-raw-to-proj-data") + +add_executable(${EXE_NAME} xreg_convert_rad_raw_to_proj_data_main.cpp) + +target_link_libraries(${EXE_NAME} PUBLIC ${XREG_EXE_LIBS_TO_LINK}) + +install(TARGETS ${EXE_NAME}) + diff --git a/apps/image_io/convert_rad_raw_to_proj_data/xreg_convert_rad_raw_to_proj_data_main.cpp b/apps/image_io/convert_rad_raw_to_proj_data/xreg_convert_rad_raw_to_proj_data_main.cpp new file mode 100644 index 0000000..d5357f3 --- /dev/null +++ b/apps/image_io/convert_rad_raw_to_proj_data/xreg_convert_rad_raw_to_proj_data_main.cpp @@ -0,0 +1,91 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregProgOptUtils.h" +#include "xregH5ProjDataIO.h" +#include "xregRadRawProj.h" + +using namespace xreg; + +int main(int argc, char* argv[]) +{ + constexpr int kEXIT_VAL_SUCCESS = 0; + constexpr int kEXIT_VAL_BAD_USE = 1; + + ProgOpts po; + + xregPROG_OPTS_SET_COMPILE_DATE(po); + + po.set_help("Converts a collection of 2D projection images from .rad/.raw format " + "into an HDF5 projection data file for use by xReg. " + "This may be used to convert projections from the Ljubljana 2D/3D datasets " + "located at: http://lit.fe.uni-lj.si/tools.php?lang=eng."); + po.set_arg_usage(" <.rad path #1> [<.rad path #2> [... <.rad path #N>]]"); + po.set_min_num_pos_args(2); + + try + { + po.parse(argc, argv); + } + catch (const ProgOpts::Exception& e) + { + std::cerr << "Error parsing command line arguments: " << e.what() << std::endl; + po.print_usage(std::cerr); + return kEXIT_VAL_BAD_USE; + } + + if (po.help_set()) + { + po.print_usage(std::cout); + po.print_help(std::cout); + return kEXIT_VAL_SUCCESS; + } + + std::ostream& vout = po.vout(); + + const std::string& dst_pd_path = po.pos_args()[0]; + + const size_type num_rad_files = po.pos_args().size() - 1; + + vout << "number of .rad files to read: " << num_rad_files << std::endl; + + ProjDataF32List pd; + pd.reserve(num_rad_files); + + for (size_type i = 0; i < num_rad_files; ++i) + { + const std::string& cur_rad_path = po.pos_args()[i+1]; + + vout << "reading .rad file #" << (i + 1) << ": " << cur_rad_path << std::endl; + + pd.push_back(ReadRawProjAsProjData(cur_rad_path)); + } + + vout << "saving to proj data HDF5..." << std::endl; + WriteProjDataH5ToDisk(pd, dst_pd_path); + + vout << "exiting..." << std::endl; + return kEXIT_VAL_SUCCESS; +} + diff --git a/apps/image_io/convert_radiograph_dicom_to_proj_data/CMakeLists.txt b/apps/image_io/convert_radiograph_dicom_to_proj_data/CMakeLists.txt new file mode 100644 index 0000000..3d1fbcc --- /dev/null +++ b/apps/image_io/convert_radiograph_dicom_to_proj_data/CMakeLists.txt @@ -0,0 +1,30 @@ +# MIT License +# +# Copyright (c) 2021 Robert Grupp +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +set(EXE_NAME "${XREG_EXE_PREFIX}convert-radiograph-dcm-to-proj-data") + +add_executable(${EXE_NAME} xreg_convert_radiograph_dcm_to_proj_data_main.cpp) + +target_link_libraries(${EXE_NAME} PUBLIC ${XREG_EXE_LIBS_TO_LINK}) + +install(TARGETS ${EXE_NAME}) + diff --git a/apps/image_io/convert_radiograph_dicom_to_proj_data/xreg_convert_radiograph_dcm_to_proj_data_main.cpp b/apps/image_io/convert_radiograph_dicom_to_proj_data/xreg_convert_radiograph_dcm_to_proj_data_main.cpp new file mode 100644 index 0000000..f7747cd --- /dev/null +++ b/apps/image_io/convert_radiograph_dicom_to_proj_data/xreg_convert_radiograph_dcm_to_proj_data_main.cpp @@ -0,0 +1,193 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregProgOptUtils.h" +#include "xregH5ProjDataIO.h" +#include "xregReadProjDataFromDICOM.h" +#include "xregStringUtils.h" + +int main(int argc, char* argv[]) +{ + using namespace xreg; + + constexpr int kEXIT_VAL_SUCCESS = 0; + constexpr int kEXIT_VAL_BAD_USE = 1; + + ProgOpts po; + + xregPROG_OPTS_SET_COMPILE_DATE(po); + + po.set_help("Convert a DICOM radiograph file into a xReg HDF5 proj. data file. " + "Support for video fluoroscopy is provided using the number of frames DICOM tag."); + po.set_arg_usage(" []"); + po.set_min_num_pos_args(2); + + po.add("src-to-det", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_DOUBLE, "src-to-det", + "Source to detector (mm) value to use ONLY when the corresponding DICOM field is not populated.") + << 1000.0; + + po.add("spacing", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_DOUBLE, "spacing", + "Image row and column spacing (mm / pixel) value to use ONLY when a suitable value may " + "not be obtained from the DICOM metadata.") + << 1.0; + + po.add("no-guess-spacing", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-guess-spacing", + "Do NOT guess pixel spacings based on other metadata values, such FOV shape and size. " + "A guess will override any value set by \"spacing\" unless the metadata needed to make " + "a guess is not available.") + << false; + + po.add("proj-frame", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_STRING, "proj-frame", + "Orientation of the projective coordinate frame to be used unless overriden by another " + "superceding flag. " + "Origin is always at the X-ray source. " + "The X axis runs parallel with the direction along an image row and is oriented in the " + "increasing column direction. " + "The Y axis runs parallel with the direction along an image column and is oriented " + "in the increasing row direction. " + "The Z axis is orthogonal to the detector plane and this flag determines the direction " + "(either towards the X-ray source or away). " + "Values: " + "\"det-neg-z\" --> Moving from the X-ray source to the detector is a movement in the " + "negative Z direction. " + "\"det-pos-z\" --> Moving from the X-ray source to the detector is a movement in the " + "positive Z direction. " + "\"auto\" --> Automatically select the orientation based on the modality field. " + "\"XA\" and \"RF\" yield \"det-neg-z\" while \"CR\" and \"DX\" yield \"det-pos-z\"") + << "auto"; + + po.add("no-proc", 'n', ProgOpts::kSTORE_TRUE, "no-proc", + "Do not perform any pre-processing to the image pixels - e.g. do NOT flip " + "or rotate the image using the DICOM FOV Rotation or FOV Horizontal Flip fields.") + << false; + + po.add("no-bayview-check", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-bayview-check", + "Do NOT inspect metadata fields to determine that the input file was created using the " + "Siemens CIOS Fusion C-arm in the JHMI Bayview lab. When this check IS performed and a " + "file is determined to have been created using the Bayview C-arm, then the extrinsic " + "transformation of the projection is populated from a previously computed transformation " + "which has an X-axis and center of rotation corresponding to the oribital rotation of " + "the C-arm.") + << false; + + po.add("pixel-type", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_STRING, "pixel-type", + "Pixel type used when saving the output projection data. Valid values are: " + "\"float\" for 32-bit floats and \"uint16\" for unsigned 16-bit integers.") + << "float"; + + po.add("fcsv-spacing", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_DOUBLE, "fcsv-spacing", + "Default (isotopic) pixel spacing to assume when parsing FCSV landmarks when no 2D pixel " + "spacing is provided by the 2D image metadata.") + << 1.0; + + try + { + po.parse(argc, argv); + } + catch (const ProgOpts::Exception& e) + { + std::cerr << "Error parsing command line arguments: " << e.what() << std::endl; + po.print_usage(std::cerr); + return kEXIT_VAL_BAD_USE; + } + + if (po.help_set()) + { + po.print_usage(std::cout); + po.print_help(std::cout); + return kEXIT_VAL_SUCCESS; + } + + std::ostream& vout = po.vout(); + + ReadProjDataFromDICOMParams read_dcm_params; + + read_dcm_params.vout = &vout; + + read_dcm_params.src_to_det_default = po.get("src-to-det"); + read_dcm_params.spacing_default = po.get("spacing"); + + read_dcm_params.guess_spacing = !po.get("no-guess-spacing"); + + read_dcm_params.fcsv_spacing_default = po.get("fcsv-spacing"); + + const std::string proj_frame_str = ToLowerCase(po.get("proj-frame").as_string()); + + if (proj_frame_str == "auto") + { + vout << "automatically selecting proj. frame Z direction using modality..." << std::endl; + } + else if (proj_frame_str == "det-neg-z") + { + vout << "forcing proj. frame Z direction to det-neg-z" << std::endl; + read_dcm_params.proj_frame = CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z; + } + else if (proj_frame_str == "det-pos-z") + { + vout << "forcing proj. frame Z direction to det-pos-z" << std::endl; + read_dcm_params.proj_frame = CameraModel::kORIGIN_AT_FOCAL_PT_DET_POS_Z; + } + else + { + std::cerr << "ERROR: Unsupported \"proj-frame\" value: " << proj_frame_str << std::endl; + return kEXIT_VAL_BAD_USE; + } + + read_dcm_params.no_proc = po.get("no-proc"); + + read_dcm_params.no_bayview_check = po.get("no-bayview-check"); + + const std::string pixel_type_str = po.get("pixel-type"); + + const std::string& src_dcm_path = po.pos_args()[0]; + const std::string& dst_pd_path = po.pos_args()[1]; + + const std::string fcsv_path = (po.pos_args().size() > 2) ? po.pos_args()[2] : std::string(); + + if (pixel_type_str == "float") + { + vout << "populating float32 proj. data from DICOM..." << std::endl; + const auto pd = ReadProjDataFromDICOMF32(src_dcm_path, fcsv_path, read_dcm_params); + + vout << " saving to proj data HDF5..." << std::endl; + WriteProjDataH5ToDisk(pd, dst_pd_path); + } + else if (pixel_type_str == "uint16") + { + vout << "populating uint16 proj. data from DICOM..." << std::endl; + const auto pd = ReadProjDataFromDICOMU16(src_dcm_path, fcsv_path, read_dcm_params); + + vout << " saving to proj data HDF5..." << std::endl; + WriteProjDataH5ToDisk(pd, dst_pd_path); + } + else + { + std::cerr << "ERROR: unsupported output pixel type: " << pixel_type_str << std::endl; + + return kEXIT_VAL_BAD_USE; + } + + return kEXIT_VAL_SUCCESS; +} + diff --git a/apps/image_io/convert_resample_dicom/xreg_convert_resample_dicom_main.cpp b/apps/image_io/convert_resample_dicom/xreg_convert_resample_dicom_main.cpp index 015d9f7..2df65b5 100644 --- a/apps/image_io/convert_resample_dicom/xreg_convert_resample_dicom_main.cpp +++ b/apps/image_io/convert_resample_dicom/xreg_convert_resample_dicom_main.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -73,8 +73,8 @@ void ReadDICOMPixelsResampleAndWriteVolume(const DICOMFIleBasicFieldsList& dcm_i using SliceType = itk::Image; using SlicePointer = typename SliceType::Pointer; - if (dcm_infos[0].sec_cap_dev_software_versions_valid && - (dcm_infos[0].sec_cap_dev_software_versions.find("Amira") != std::string::npos)) + if (dcm_infos[0].sec_cap_dev_software_versions && + (dcm_infos[0].sec_cap_dev_software_versions->find("Amira") != std::string::npos)) { std::cerr << "WARNING: THIS SERIES WAS EXPORTED FROM AMIRA - " "LPS DIRECTIONS MAY NOT BE CORRECT IF RESAMPLING/WARPING WAS " @@ -95,98 +95,109 @@ void ReadDICOMPixelsResampleAndWriteVolume(const DICOMFIleBasicFieldsList& dcm_i Pt3 out_of_plane_dir = in_plane_col_dir.cross(in_plane_row_dir); { - const std::string pat_pos_str = StringStrip(dcm_infos[0].pat_pos); - const std::string pat_orient_x = StringStrip(dcm_infos[0].pat_orient[0]); - const std::string pat_orient_y = StringStrip(dcm_infos[0].pat_orient[1]); + const bool has_pat_orient_strs = dcm_infos[0].pat_orient.has_value(); - if (pat_pos_str == "HFS") // head first supine, most common + const std::string pat_orient_x = has_pat_orient_strs ? + StringStrip((*dcm_infos[0].pat_orient)[0]) : std::string(); + const std::string pat_orient_y = has_pat_orient_strs ? + StringStrip((*dcm_infos[0].pat_orient)[1]) : std::string(); + + if (dcm_infos[0].pat_pos) { - // This should have pixels already in LPS - if ((pat_orient_x != "L") || (pat_orient_y != "P")) + const std::string pat_pos_str = StringStrip(*dcm_infos[0].pat_pos); + + if (pat_pos_str == "HFS") // head first supine, most common { - std::cerr << "WARNING: Patient position (" << pat_pos_str - << ") is inconsistent with patient orientation (" - << pat_orient_x << " , " << pat_orient_y << ')' << std::endl; + // This should have pixels already in LPS + if (has_pat_orient_strs && ((pat_orient_x != "L") || (pat_orient_y != "P"))) + { + std::cerr << "WARNING: Patient position (" << pat_pos_str + << ") is inconsistent with patient orientation (" + << pat_orient_x << " , " << pat_orient_y << ')' << std::endl; + } } - } - else if(pat_pos_str == "HFP") // head first prone, have seen a few of these - { - // This should have pixels in RAS - if ((pat_orient_x != "R") || (pat_orient_y != "A")) + else if(pat_pos_str == "HFP") // head first prone, have seen a few of these { - std::cerr << "WARNING: Patient position (" << pat_pos_str - << ") is inconsistent with patient orientation (" - << pat_orient_x << " , " << pat_orient_y << ')' << std::endl; + // This should have pixels in RAS + if (has_pat_orient_strs && ((pat_orient_x != "R") || (pat_orient_y != "A"))) + { + std::cerr << "WARNING: Patient position (" << pat_pos_str + << ") is inconsistent with patient orientation (" + << pat_orient_x << " , " << pat_orient_y << ')' << std::endl; + } + } + //else if (pat_pos_str == "FFS") // Feet first supine + else + { + std::cerr << "WARNING: UNCOMMON PATIENT POSITON STRING: " + << pat_pos_str << std::endl; } - } - //else if (pat_pos_str == "FFS") // Feet first supine - else - { - std::cerr << "WARNING: UNCOMMON PATIENT POSITON STRING: " - << pat_pos_str << std::endl; } - Pt3 std_x_axis; - std_x_axis[0] = 1; - std_x_axis[1] = 0; - std_x_axis[2] = 0; + if (has_pat_orient_strs) + { + Pt3 std_x_axis; + std_x_axis[0] = 1; + std_x_axis[1] = 0; + std_x_axis[2] = 0; - Pt3 std_y_axis; - std_y_axis[0] = 0; - std_y_axis[1] = 1; - std_y_axis[2] = 0; + Pt3 std_y_axis; + std_y_axis[0] = 0; + std_y_axis[1] = 1; + std_y_axis[2] = 0; - if (pat_orient_x == "L") - { - if ((in_plane_col_dir - std_x_axis).norm() > 1.0e-6) + if (pat_orient_x == "L") + { + if ((in_plane_col_dir - std_x_axis).norm() > 1.0e-6) + { + std::cerr << "WARNING: PATIENT ORIENTATION OF X-AXIS (" << pat_orient_x + << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" + << in_plane_col_dir[0] << ',' << in_plane_col_dir[1] + << ',' << in_plane_col_dir[2] << ')' << std::endl; + } + } + else if (pat_orient_x == "R") { - std::cerr << "WARNING: PATIENT ORIENTATION OF X-AXIS (" << pat_orient_x - << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" - << in_plane_col_dir[0] << ',' << in_plane_col_dir[1] - << ',' << in_plane_col_dir[2] << ')' << std::endl; + if ((in_plane_col_dir + std_x_axis).norm() > 1.0e-6) + { + std::cerr << "WARNING: PATIENT ORIENTATION OF X-AXIS (" << pat_orient_x + << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" + << in_plane_col_dir[0] << ',' << in_plane_col_dir[1] + << ',' << in_plane_col_dir[2] << ')' << std::endl; + } } - } - else if (pat_orient_x == "R") - { - if ((in_plane_col_dir + std_x_axis).norm() > 1.0e-6) + else { - std::cerr << "WARNING: PATIENT ORIENTATION OF X-AXIS (" << pat_orient_x - << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" - << in_plane_col_dir[0] << ',' << in_plane_col_dir[1] - << ',' << in_plane_col_dir[2] << ')' << std::endl; + std::cerr << "WARNING: UNCOMMON PATIENT ORIENTATION OF X-AXIS (" + << pat_orient_x << ')' << std::endl; } - } - else - { - std::cerr << "WARNING: UNCOMMON PATIENT ORIENTATION OF X-AXIS (" - << pat_orient_x << ')' << std::endl; - } - if (pat_orient_y == "P") - { - if ((in_plane_row_dir - std_y_axis).norm() > 1.0e-6) + if (pat_orient_y == "P") { - std::cerr << "WARNING: PATIENT ORIENTATION OF Y-AXIS (" << pat_orient_y - << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" - << in_plane_row_dir[0] << ',' << in_plane_row_dir[1] - << ',' << in_plane_row_dir[2] << ')' << std::endl; + if ((in_plane_row_dir - std_y_axis).norm() > 1.0e-6) + { + std::cerr << "WARNING: PATIENT ORIENTATION OF Y-AXIS (" << pat_orient_y + << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" + << in_plane_row_dir[0] << ',' << in_plane_row_dir[1] + << ',' << in_plane_row_dir[2] << ')' << std::endl; + } } - } - else if (pat_orient_y == "A") - { - if ((in_plane_row_dir + std_y_axis).norm() > 1.0e-6) + else if (pat_orient_y == "A") + { + if ((in_plane_row_dir + std_y_axis).norm() > 1.0e-6) + { + std::cerr << "WARNING: PATIENT ORIENTATION OF Y-AXIS (" << pat_orient_y + << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" + << in_plane_row_dir[0] << ',' << in_plane_row_dir[1] + << ',' << in_plane_row_dir[2] << ')' << std::endl; + } + } + else { - std::cerr << "WARNING: PATIENT ORIENTATION OF Y-AXIS (" << pat_orient_y - << ") IS NOT CONSISTENT WITH DIRECTION VECTOR (" - << in_plane_row_dir[0] << ',' << in_plane_row_dir[1] - << ',' << in_plane_row_dir[2] << ')' << std::endl; + std::cerr << "WARNING: UNCOMMON PATIENT ORIENTATION OF Y-AXIS (" + << pat_orient_y << ')' << std::endl; } } - else - { - std::cerr << "WARNING: UNCOMMON PATIENT ORIENTATION OF Y-AXIS (" - << pat_orient_y << ')' << std::endl; - } } if (out_of_plane_dir[out_of_plane_axis_dim] < 0) @@ -735,12 +746,23 @@ int main( int argc, char* argv[] ) "Limit processing to a specific series UID.") << ""; + po.add("no-name-w-pat-id", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-name-w-pat-id", + "Do NOT include the patient ID for naming output files. Use care with this option as " + "datasets may have different patient IDs with the same patient name, which could " + "possibily result in several series mapping to the same output file. This could result " + "in not all series being written out to disk WITHOUT WARNING.") + << false; + + po.add("name-w-pat-name", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "name-w-pat-name", + "Include patient name for naming output files.") + << false; + po.add("name-w-desc", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "name-w-desc", - "Use study and series descriptions for naming output files.") + "Include study and series descriptions for naming output files.") << false; po.add("name-w-conv", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "name-w-conv", - "Use convolution kernel for naming output files (useful for distinguishing between " + "Include convolution kernel for naming output files (useful for distinguishing between " "soft-tissue and bone reconstructions).") << false; @@ -755,6 +777,22 @@ int main( int argc, char* argv[] ) "or a sequence of 2D images.)") << false; + po.add("exclude-derived", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "exclude-derived", + "Do NOT include DERIVED images (e.g. images with pixel values derived from other images).") + << false; + + po.add("exclude-secondary", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "exclude-secondary", + "Do NOT include SECONDARY images (e.g. images created after initial exam).") + << false; + + po.add("modalities", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_STRING, "modalities", + "A comma delimited string indicating which modalities should be considered. " + "An empty string (the default) considers all modalities. " + "This is case-sensitive - upper case strings are standard. " + "Examples: \"CT\" will only convert CT datasets, " + "\"CT,MR\" will convert CT and MR datasets.") + << ""; + po.add("pat-lut", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_STRING, "pat-lut", "Path to a CSV file which serves as a LUT between DICOM patient IDs and strings " "used to prefix output file names. The first column are the DICOM patient IDs " @@ -876,23 +914,27 @@ int main( int argc, char* argv[] ) const std::string single_series_uid = po.get("series-uid"); const bool limit_single_series = !single_series_uid.empty(); + const bool use_pat_id_for_name = !po.get("no-name-w-pat-id").as_bool(); + + const bool use_pat_name_for_name = po.get("name-w-pat-name"); + const bool use_desc_for_name = po.get("name-w-desc"); const bool use_conv_for_name = po.get("name-w-conv"); - if ((int(use_desc_for_name) + int(use_conv_for_name)) > 1) - { - std::cerr << "ERROR: can only use one flag for choosing DICOM fields for output naming!" - << " Choose one of \"name-w-desc\" or \"name-w-conv\"" << std::endl; - - return kEXIT_VAL_BAD_USE; - } - const bool inc_localizers = po.get("include-localizer"); const bool inc_multiframe_files = po.get("include-multiframe-files"); + const bool exclude_derived = po.get("exclude-derived"); + + const bool exclude_secondary = po.get("exclude-secondary"); + const std::string pat_id_lut_path = po.get("pat-lut"); + const std::string modalities_to_consider_str = po.get("modalities"); + + const auto modalities_to_consider = StringSplit(modalities_to_consider_str, ","); + const bool use_pat_id_lut = !pat_id_lut_path.empty(); vout << "Inputs: " @@ -933,10 +975,15 @@ int main( int argc, char* argv[] ) << "\n Single Patient ID: " << (limit_single_pat ? single_pat_id.c_str() : "(All Patients)") << "\n Single Study UID: " << (limit_single_study ? single_study_uid.c_str() : "(All Studies)") << "\n Single Series UID: " << (limit_single_series ? single_series_uid.c_str() : "(All Series)") + << "\n Use Patient ID for Name: " << use_pat_id_for_name + << "\n Use Patient Name for Name: " << use_pat_name_for_name << "\n Use Study/Series Desc. Name: " << use_desc_for_name << "\n Use Conv. Kern. Name: " << use_conv_for_name << "\n Include Localizers: " << inc_localizers << "\n Include Multi-Frame Files: " << inc_multiframe_files + << "\n Exclude Derived Files: " << exclude_derived + << "\n Exclude Secondary Files: " << exclude_secondary + << "\n Modalities Considered: " << (modalities_to_consider_str.empty() ? "" : modalities_to_consider_str.c_str()) << "\n Use Patient ID LUT: " << use_pat_id_lut << "\n-----------------------------------------------------\n" << std::endl; @@ -965,7 +1012,9 @@ int main( int argc, char* argv[] ) { vout << "reading directory tree and organizing..." << std::endl; GetOrgainizedDICOMInfos(input_root_dir, &org_dcm, - inc_localizers, inc_multiframe_files); + inc_localizers, inc_multiframe_files, + !exclude_secondary, !exclude_derived, + modalities_to_consider); } else { @@ -989,8 +1038,8 @@ int main( int argc, char* argv[] ) org_dcm.root_dir = input_root_dir; // manually assign all DICOM files to a single patient and study - DICOMFIleBasicFields dcm_info; - ReadDICOMFileBasicFields(dcm_str_paths[0], &dcm_info); + DICOMFIleBasicFields dcm_info = ReadDICOMFileBasicFields(dcm_str_paths[0]); + org_dcm.patient_infos[dcm_info.patient_id] [dcm_info.study_uid][dcm_info.series_uid] = dcm_str_paths; } @@ -1070,13 +1119,16 @@ int main( int argc, char* argv[] ) vout << std::endl; } - auto get_valid_name_for_path = [file_name_char_remap] (const bool valid, const std::string& s) + // Replace invalid/undesired characters in a path string using the file_name_char_remap LUT. + // When the optional variable for the input string indicates the value is not available, then + // and empty string is returned + auto get_valid_name_for_path = [file_name_char_remap] (const boost::optional& s) { std::string ss; - if (valid) + if (s) { - ss = MapChars(s, file_name_char_remap); + ss = MapChars(*s, file_name_char_remap); } return ss; @@ -1158,7 +1210,7 @@ int main( int argc, char* argv[] ) for (size_type dcm_idx = 0; dcm_idx < num_dcm; ++dcm_idx) { - ReadDICOMFileBasicFields(paths[dcm_idx], &dcm_infos[dcm_idx]); + dcm_infos[dcm_idx] = ReadDICOMFileBasicFields(paths[dcm_idx]); } ReorderAndCheckDICOMInfos dcm_check_reorder; @@ -1241,31 +1293,38 @@ int main( int argc, char* argv[] ) if (!one_image_input) { - std::string dst_file_name; - - const std::string pat_id_for_filename = get_valid_name_for_path(true, pat_id_str_remap); - - if (use_desc_for_name) + std::stringstream dst_file_name_ss; + + if (use_pat_id_for_name) { - dst_file_name = fmt::format("{}_{:02d}_{}_{}.{}", - pat_id_for_filename, cur_img_idx, - get_valid_name_for_path(first_dcm.study_desc_valid, first_dcm.study_desc), - get_valid_name_for_path(first_dcm.series_desc_valid, first_dcm.series_desc), - out_file_ext); + dst_file_name_ss << get_valid_name_for_path(pat_id_str_remap); + dst_file_name_ss << '_'; } - else if (use_conv_for_name) + + if (use_pat_name_for_name) + { + dst_file_name_ss << get_valid_name_for_path(first_dcm.patient_name) + << '_'; + } + + dst_file_name_ss << fmt::format("{:02d}", cur_img_idx); + + if (use_desc_for_name) { - dst_file_name = fmt::format("{}_{:02d}_{}.{}", - pat_id_for_filename, cur_img_idx, - get_valid_name_for_path(first_dcm.conv_kernel_valid, first_dcm.conv_kernel), - out_file_ext); + dst_file_name_ss << '_' + << get_valid_name_for_path(first_dcm.study_desc) + << get_valid_name_for_path(first_dcm.series_desc); } - else + + if (use_conv_for_name) { - dst_file_name = fmt::sprintf("%s_%02lu.%s", pat_id_str_remap, cur_img_idx, out_file_ext); + dst_file_name_ss << '_' + << get_valid_name_for_path(first_dcm.conv_kernel); } + + dst_file_name_ss << '.' << out_file_ext; - dst_file_path += dst_file_name; + dst_file_path += dst_file_name_ss.str(); } else if (dst_file_path.file_extension().empty()) { diff --git a/apps/image_io/convert_sta_raw_to_itk/CMakeLists.txt b/apps/image_io/convert_sta_raw_to_itk/CMakeLists.txt new file mode 100644 index 0000000..bad095a --- /dev/null +++ b/apps/image_io/convert_sta_raw_to_itk/CMakeLists.txt @@ -0,0 +1,30 @@ +# MIT License +# +# Copyright (c) 2021 Robert Grupp +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +set(EXE_NAME "${XREG_EXE_PREFIX}convert-sta-raw-to-itk") + +add_executable(${EXE_NAME} xreg_convert_sta_raw_to_itk_main.cpp) + +target_link_libraries(${EXE_NAME} PUBLIC ${XREG_EXE_LIBS_TO_LINK}) + +install(TARGETS ${EXE_NAME}) + diff --git a/apps/image_io/convert_sta_raw_to_itk/xreg_convert_sta_raw_to_itk_main.cpp b/apps/image_io/convert_sta_raw_to_itk/xreg_convert_sta_raw_to_itk_main.cpp new file mode 100644 index 0000000..a4d5e69 --- /dev/null +++ b/apps/image_io/convert_sta_raw_to_itk/xreg_convert_sta_raw_to_itk_main.cpp @@ -0,0 +1,79 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregProgOptUtils.h" +#include "xregITKIOUtils.h" +#include "xregStaRawVol.h" + +int main(int argc, char* argv[]) +{ + using namespace xreg; + + constexpr int kEXIT_VAL_SUCCESS = 0; + constexpr int kEXIT_VAL_BAD_USE = 1; + + ProgOpts po; + + xregPROG_OPTS_SET_COMPILE_DATE(po); + + po.set_help("Convert a .sta/.raw volume to an ITK compatible volume file (e.g. .nii.gz). " + "The .sta/.raw format is used in the Ljubljana 2D/3D datasets " + "(http://lit.fe.uni-lj.si/tools.php?lang=eng)."); + po.set_arg_usage(" "); + + po.set_min_num_pos_args(2); + + try + { + po.parse(argc, argv); + } + catch (const xreg::ProgOpts::Exception& e) + { + std::cerr << "Error parsing command line arguments: " << e.what() << std::endl; + po.print_usage(std::cerr); + return kEXIT_VAL_BAD_USE; + } + + if (po.help_set()) + { + po.print_usage(std::cout); + po.print_help(std::cout); + return kEXIT_VAL_SUCCESS; + } + + const std::string input_sta_path = po.pos_args()[0]; + const std::string output_vol_path = po.pos_args()[1]; + + std::ostream& vout = po.vout(); + + vout << "reading .sta/.raw data..." << std::endl; + auto vol = ReadStaRawVol(input_sta_path); + + vout << "writing to output ITK volume file..." << std::endl; + WriteITKImageToDisk(vol.GetPointer(), output_vol_path); + + vout << "exiting..." << std::endl; + + return kEXIT_VAL_SUCCESS; +} diff --git a/apps/image_io/extract_nii_from_proj_data/xreg_extract_nii_from_proj_data_main.cpp b/apps/image_io/extract_nii_from_proj_data/xreg_extract_nii_from_proj_data_main.cpp index c9e67aa..ab6b76c 100644 --- a/apps/image_io/extract_nii_from_proj_data/xreg_extract_nii_from_proj_data_main.cpp +++ b/apps/image_io/extract_nii_from_proj_data/xreg_extract_nii_from_proj_data_main.cpp @@ -38,6 +38,7 @@ void ProcessAndSave(itk::Image* img, const std::string& path, const boost::optional& rot_to_pat_up, const double ds_factor, + const bool discard_geom, std::ostream& vout) { using Img = itk::Image; @@ -60,6 +61,18 @@ void ProcessAndSave(itk::Image* img, tmp_img = DownsampleImage(img_to_save, ds_factor); img_to_save = tmp_img.GetPointer(); } + + if (discard_geom) + { + const std::array spacing = { 1.0, 1.0 }; + const std::array origin = { 0.0, 0.0 }; + typename Img::DirectionType dir_mat; + dir_mat.SetIdentity(); + + img_to_save->SetSpacing(spacing.data()); + img_to_save->SetOrigin(origin.data()); + img_to_save->SetDirection(dir_mat); + } vout << " writing image..." << std::endl; WriteITKImageToDisk(img_to_save, path); @@ -96,6 +109,12 @@ int main(int argc, char* argv[]) po.add("ds-factor", 'd', ProgOpts::kSTORE_DOUBLE, "ds-factor", "Downsampling factor applied to projection data") << 1.0; + + po.add("no-geom", ProgOpts::kNO_SHORT_FLAG, ProgOpts::kSTORE_TRUE, "no-geom", + "Discard any geometric metadata (e.g. pixel spacings, origin, orientation) when writing " + "NIFTI files. The output file will have pixel spacings equal to one, identity orientation " + "and zero origin.") + << false; try { @@ -121,6 +140,8 @@ int main(int argc, char* argv[]) const double proj_ds_factor = po.get("ds-factor"); + const bool discard_geom = po.get("no-geom"); + const std::string proj_data_path = po.pos_args()[0]; const std::string nii_prefix = po.pos_args()[1]; const std::string proj_inds_str = (po.pos_args().size()) > 2 ? po.pos_args()[2] : std::string(); @@ -191,19 +212,19 @@ int main(int argc, char* argv[]) { auto img = pd_reader.read_proj_F32(src_proj_idx); - ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, vout); + ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, discard_geom, vout); } else if (scalar_type == kPROJ_DATA_TYPE_UINT16) { auto img = pd_reader.read_proj_U16(src_proj_idx); - ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, vout); + ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, discard_geom, vout); } else if (scalar_type == kPROJ_DATA_TYPE_UINT8) { auto img = pd_reader.read_proj_U8(src_proj_idx); - ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, vout); + ProcessAndSave(img.GetPointer(), dst_path, rot_to_pat_up, proj_ds_factor, discard_geom, vout); } } } diff --git a/apps/image_io/make_video_from_image_dir/CMakeLists.txt b/apps/image_io/make_video_from_image_dir/CMakeLists.txt new file mode 100644 index 0000000..17c17c5 --- /dev/null +++ b/apps/image_io/make_video_from_image_dir/CMakeLists.txt @@ -0,0 +1,30 @@ +# MIT License +# +# Copyright (c) 2021 Robert Grupp +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +set(EXE_NAME "${XREG_EXE_PREFIX}vid-from-img-dir") + +add_executable(${EXE_NAME} xreg_make_video_from_image_dir_main.cpp) + +target_link_libraries(${EXE_NAME} PUBLIC ${XREG_EXE_LIBS_TO_LINK}) + +install(TARGETS ${EXE_NAME}) + diff --git a/apps/image_io/make_video_from_image_dir/xreg_make_video_from_image_dir_main.cpp b/apps/image_io/make_video_from_image_dir/xreg_make_video_from_image_dir_main.cpp new file mode 100644 index 0000000..c97d726 --- /dev/null +++ b/apps/image_io/make_video_from_image_dir/xreg_make_video_from_image_dir_main.cpp @@ -0,0 +1,116 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregCommon.h" +#include "xregProgOptUtils.h" +#include "xregWriteVideo.h" + +int main(int argc, char* argv[]) +{ + using namespace xreg; + + constexpr int kEXIT_VAL_SUCCESS = 0; + constexpr int kEXIT_VAL_BAD_USE = 1; + + ProgOpts po; + + xregPROG_OPTS_SET_COMPILE_DATE(po); + + po.set_help("Create a video from a directory of image files, which are used as frames."); + po.set_arg_usage(" [ [ ... []]]"); + po.set_min_num_pos_args(2); + + po.add("sort", 's', ProgOpts::kSTORE_TRUE, "sort", "Sort frame order by lexographic file path order.") + << false; + + po.add("fps", 'f', ProgOpts::kSTORE_DOUBLE, "fps", "The frame rate of the output movie.") + << 1.0; + + po.add("len", 'l', ProgOpts::kSTORE_DOUBLE, "len", + "The length of the output movie in seconds; this argument will " + "take precedent over the fps argument."); + + try + { + po.parse(argc, argv); + } + catch (const ProgOpts::Exception& e) + { + std::cerr << "Error parsing command line arguments: " << e.what() << std::endl; + po.print_usage(std::cerr); + return kEXIT_VAL_BAD_USE; + } + + if (po.help_set()) + { + po.print_usage(std::cout); + po.print_help(std::cout); + return kEXIT_VAL_SUCCESS; + } + + std::ostream& vout = po.vout(); + + const size_type num_cmd_args = po.pos_args().size(); + + const size_type num_exts_on_cmd_line = num_cmd_args - 2; + + const std::string src_dir = po.pos_args()[0]; + const std::string dst_mov = po.pos_args()[1]; + + double fps_or_len = po.get("fps"); + + const bool sort_paths = po.get("sort"); + + const bool video_len_passed = po.has("len"); + + double video_len_secs = video_len_passed ? po.get("len").as_double() : 0.0; + + vout << num_exts_on_cmd_line << " extensions passed on the command line..." << std::endl; + + std::vector img_exts; + + if (num_exts_on_cmd_line) + { + img_exts.insert(img_exts.end(), po.pos_args().begin() + 2, po.pos_args().end()); + } + else + { + vout << "using default image extension of .png" << std::endl; + + img_exts = { ".png" }; + } + + if (video_len_passed) + { + vout << "overriding FPS with video length argument..." << std::endl; + fps_or_len = po.get("len"); + } + + vout << "processing directory contents and creating video..." << std::endl; + WriteDirOfImagesToVideo(dst_mov, src_dir, sort_paths, img_exts, fps_or_len, + !video_len_passed); + + return kEXIT_VAL_SUCCESS; +} + diff --git a/apps/image_io/regi2d3d_replay/xreg_regi2d3d_replay_main.cpp b/apps/image_io/regi2d3d_replay/xreg_regi2d3d_replay_main.cpp index a78df84..5c5681d 100644 --- a/apps/image_io/regi2d3d_replay/xreg_regi2d3d_replay_main.cpp +++ b/apps/image_io/regi2d3d_replay/xreg_regi2d3d_replay_main.cpp @@ -331,7 +331,7 @@ int main(int argc, char* argv[]) // check to see if a subset of views are used... size_type num_views_used = 0; - // this is not ideal, but it gets us by for know + // this is not ideal, but it gets us by for now // ensure that if a subset of views is specified, then the same views are used in all registrations { auto projs_used_are_same = [] (const IndexList& i1, const IndexList& i2) @@ -484,7 +484,7 @@ int main(int argc, char* argv[]) if (std::abs(1.0 - ds_factor) > 1.0e-6) { vout << "downsampling fixed images..." << std::endl; - fixed_proj_data = DownsampleProjData(fixed_proj_data, ds_factor); + fixed_proj_data = DownsampleProjData(fixed_proj_data, ds_factor, true); } else { @@ -1023,7 +1023,7 @@ int main(int argc, char* argv[]) frame_file_names.reserve(tot_num_projs); frame_titles.reserve(tot_num_projs); - const int tile_border_thickness = std::min(1, static_cast(std::lround(10 * ds_factor))); + const int tile_border_thickness = std::max(1, static_cast(std::lround(10 * ds_factor))); for (size_type lvl = 0; lvl < num_levels; ++lvl) { diff --git a/apps/image_io/remap_and_tile_proj_data/xreg_remap_tile_proj_data_main.cpp b/apps/image_io/remap_and_tile_proj_data/xreg_remap_tile_proj_data_main.cpp index cef8f75..30390af 100644 --- a/apps/image_io/remap_and_tile_proj_data/xreg_remap_tile_proj_data_main.cpp +++ b/apps/image_io/remap_and_tile_proj_data/xreg_remap_tile_proj_data_main.cpp @@ -299,7 +299,10 @@ int main(int argc, char* argv[]) std::lround(3 * (((((proj_ds_factor < 1) ? -1 : 1) * 0.75) * (proj_ds_factor - 1) * (proj_ds_factor - 1)) + 1))); - border_thickness = std::max(1l, std::lround(proj_ds_factor * border_thickness)); + if (border_thickness > 0) + { + border_thickness = std::max(1l, std::lround(proj_ds_factor * border_thickness)); + } } const std::string proj_data_path = po.pos_args()[0]; diff --git a/apps/image_io/remap_dicom_dir/CMakeLists.txt b/apps/image_io/remap_dicom_dir/CMakeLists.txt new file mode 100644 index 0000000..5db892d --- /dev/null +++ b/apps/image_io/remap_dicom_dir/CMakeLists.txt @@ -0,0 +1,30 @@ +# MIT License +# +# Copyright (c) 2021 Robert Grupp +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +set(EXE_NAME "${XREG_EXE_PREFIX}remap-dcm-dir") + +add_executable(${EXE_NAME} xreg_remap_dcm_dir_main.cpp) + +target_link_libraries(${EXE_NAME} PUBLIC ${XREG_EXE_LIBS_TO_LINK}) + +install(TARGETS ${EXE_NAME}) + diff --git a/apps/image_io/remap_dicom_dir/xreg_remap_dcm_dir_main.cpp b/apps/image_io/remap_dicom_dir/xreg_remap_dcm_dir_main.cpp new file mode 100644 index 0000000..dd57c88 --- /dev/null +++ b/apps/image_io/remap_dicom_dir/xreg_remap_dcm_dir_main.cpp @@ -0,0 +1,172 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include + +#include "xregFilesystemUtils.h" +#include "xregProgOptUtils.h" +#include "xregReadProjDataFromDICOM.h" +#include "xregITKIOUtils.h" + +int main(int argc, char* argv[]) +{ + using namespace xreg; + + constexpr int kEXIT_VAL_SUCCESS = 0; + constexpr int kEXIT_VAL_BAD_USE = 1; + + ProgOpts po; + + xregPROG_OPTS_SET_COMPILE_DATE(po); + + po.set_help("Given a directory of DICOM files, writes remapped versions (for display) of each " + "file to an output directory. If the output directory does not exist, it is created. " + "For DICOM files with more than one frame, each frame is saved as a separate remapped " + "file with a \"_\" string appended before the extension (For example " + "multi_frame_dicom.dcm --> { multi_frame_dicom_0.png, ..., multi_frame_dicom_N.png })."); + po.set_arg_usage(" "); + po.set_min_num_pos_args(2); + + po.add("no-proc", 'n', ProgOpts::kSTORE_TRUE, "no-proc", + "Do not perform any pre-processing to the image pixels - e.g. do NOT flip " + "or rotate the image using the DICOM FOV Rotation or FOV Horizontal Flip fields.") + << false; + + po.add("ds", 'd', ProgOpts::kSTORE_DOUBLE, "ds", + "Downsample factor in each 2D dimension. " + "0.25 --> 2x downsampling in each dimension. " + "1 --> no downsampling. " + "2 --> 2x upsampling in each dimension.") + << 1.0; + + po.add("ext", 'e', ProgOpts::kSTORE_STRING, "ext", + "Extension (and file format) to use for remapped images.") + << "png"; + + try + { + po.parse(argc, argv); + } + catch (const ProgOpts::Exception& e) + { + std::cerr << "Error parsing command line arguments: " << e.what() << std::endl; + po.print_usage(std::cerr); + return kEXIT_VAL_BAD_USE; + } + + if (po.help_set()) + { + po.print_usage(std::cout); + po.print_help(std::cout); + return kEXIT_VAL_SUCCESS; + } + + std::ostream& vout = po.vout(); + + ReadProjDataFromDICOMParams read_dcm_params; + + read_dcm_params.vout = &vout; + + read_dcm_params.no_proc = po.get("no-proc"); + + const double ds_factor = po.get("ds"); + + const bool do_ds = std::abs(1.0 - ds_factor) > 0.001; + + const std::string out_ext = po.get("ext"); + + const std::string& src_dcm_dir_path = po.pos_args()[0]; + const std::string& dst_remap_dir_path = po.pos_args()[1]; + + Path dst_dir_path_obj(dst_remap_dir_path); + + if (!dst_dir_path_obj.exists()) + { + vout << "creating output directory..." << std::endl; + MakeDirRecursive(dst_remap_dir_path); + } + else if (!dst_dir_path_obj.is_dir()) + { + std::cerr << "ERROR: Destination directory is NOT a directory!" << std::endl; + return kEXIT_VAL_BAD_USE; + } + + PathList dcm_paths; + Path(src_dcm_dir_path).get_dir_contents(&dcm_paths); + vout << "number of entries in source dir: " << dcm_paths.size() << std::endl; + + using PDList = decltype(ReadProjDataFromDICOMF32("")); + + for (const auto& cur_path : dcm_paths) + { + PDList pd; + + try + { + pd = ReadProjDataFromDICOMF32(cur_path.string(), read_dcm_params); + } + catch (...) + { + vout << "error reading... assuming is not dicom..." << std::endl; + } + + if (!pd.empty()) + { + if (do_ds) + { + pd = DownsampleProjData(pd, ds_factor); + } + + const std::string src_filename_wo_ext = std::get<0>(cur_path.filename().split_ext()); + + const size_type num_frames = pd.size(); + + if (num_frames == 1) + { + WriteITKImageToDisk(ITKImageRemap8bpp(pd[0].img.GetPointer()).GetPointer(), + fmt::format("{}/{}.{}", dst_remap_dir_path, src_filename_wo_ext, out_ext)); + } + else + { + std::stringstream ss; + + ss << "{}/{}_{:0" + << static_cast(std::log10(num_frames) + 1) + << "d}.{}"; + + const std::string fmt_str = ss.str(); + + for (size_type frame_idx = 0; frame_idx < num_frames; ++frame_idx) + { + WriteITKImageToDisk(ITKImageRemap8bpp(pd[frame_idx].img.GetPointer()).GetPointer(), + fmt::format(fmt_str, + dst_remap_dir_path, src_filename_wo_ext, frame_idx, out_ext)); + } + } + } + } + + return kEXIT_VAL_SUCCESS; +} + diff --git a/apps/image_io/report_dicom/xreg_report_dicom_main.cpp b/apps/image_io/report_dicom/xreg_report_dicom_main.cpp index 4d7cbae..756c2cc 100644 --- a/apps/image_io/report_dicom/xreg_report_dicom_main.cpp +++ b/apps/image_io/report_dicom/xreg_report_dicom_main.cpp @@ -36,6 +36,7 @@ #include "xregProgOptUtils.h" #include "xregFilesystemUtils.h" +#include "xregDICOMUtils.h" int main( int argc, char* argv[] ) { @@ -61,6 +62,14 @@ int main( int argc, char* argv[] ) "positional arguments."); po.set_arg_usage(" [[[tag 1] tag2] ... tag N]"); + po.add("skip-full-print", 's', ProgOpts::kSTORE_TRUE, "skip-full-print", + "Skip the raw printout of all DICOM fields.") + << false; + + po.add("print-xreg-fields", 'x', ProgOpts::kSTORE_TRUE, "print-xreg-fields", + "Print basic DICOM fields as stored by the DICOMFIleBasicFields struct.") + << false; + po.set_min_num_pos_args(1); try @@ -81,88 +90,102 @@ int main( int argc, char* argv[] ) return kEXIT_VAL_SUCCESS; } + const bool skip_full_print = po.get("skip-full-print"); + + const bool print_xreg_fields = po.get("print-xreg-fields"); + using TagSetType = std::set; TagSetType in_tags(po.pos_args().begin() + 1, po.pos_args().end()); const bool limit_tags = !in_tags.empty(); std::string img_path = po.pos_args()[0]; - // First check for a directory as source input, if true then assume we are - // dealing with DICOM - itk::Directory::Pointer src_dir = itk::Directory::New(); - if (src_dir->Load(img_path.c_str())) + if (!skip_full_print) { - int good_file_idx = -1; + // First check for a directory as source input, if true then assume we are + // dealing with DICOM + itk::Directory::Pointer src_dir = itk::Directory::New(); + if (src_dir->Load(img_path.c_str())) + { + int good_file_idx = -1; - const int num_files = src_dir->GetNumberOfFiles(); + const int num_files = src_dir->GetNumberOfFiles(); - for (int file_idx = 0; file_idx < num_files; ++file_idx) - { - if (src_dir->GetFile(file_idx)[0] != '.') + for (int file_idx = 0; file_idx < num_files; ++file_idx) { - good_file_idx = file_idx; - break; + if (src_dir->GetFile(file_idx)[0] != '.') + { + good_file_idx = file_idx; + break; + } } - } - if (good_file_idx >= 0) - { - img_path += kFILE_SEPS[0]; - img_path += src_dir->GetFile(good_file_idx); - } - else - { - std::cerr << "ERROR: could not find a file to read!" << std::endl; - return kEXIT_VAL_EMPTY_DIR; + if (good_file_idx >= 0) + { + img_path += kFILE_SEPS[0]; + img_path += src_dir->GetFile(good_file_idx); + } + else + { + std::cerr << "ERROR: could not find a file to read!" << std::endl; + return kEXIT_VAL_EMPTY_DIR; + } } - } - using ImageIOType = itk::GDCMImageIO; + using ImageIOType = itk::GDCMImageIO; - ImageIOType::Pointer dcm_io = ImageIOType::New(); + ImageIOType::Pointer dcm_io = ImageIOType::New(); - using ReaderType = itk::ImageFileReader; - ReaderType::Pointer reader = ReaderType::New(); - reader->SetImageIO(dcm_io); - reader->SetFileName(img_path); + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetImageIO(dcm_io); + reader->SetFileName(img_path); - try - { - reader->Update(); - } - catch (itk::ExceptionObject& e) - { - std::cerr << "ERROR: reading file:\n" << e << std::endl; - return kEXIT_VAL_IO_FAILURE; - } - - using DictionaryType = itk::MetaDataDictionary; - - const DictionaryType& dict = dcm_io->GetMetaDataDictionary(); + try + { + reader->Update(); + } + catch (itk::ExceptionObject& e) + { + std::cerr << "ERROR: reading file:\n" << e << std::endl; + return kEXIT_VAL_IO_FAILURE; + } - DictionaryType::ConstIterator dict_end_it = dict.End(); + using DictionaryType = itk::MetaDataDictionary; - for (DictionaryType::ConstIterator dict_it = dict.Begin(); dict_it != dict_end_it; ++dict_it) - { - using MetaDataStringType = itk::MetaDataObject; + const DictionaryType& dict = dcm_io->GetMetaDataDictionary(); - MetaDataStringType::Pointer entry = dynamic_cast(dict_it->second.GetPointer()); + DictionaryType::ConstIterator dict_end_it = dict.End(); - if (entry) + for (DictionaryType::ConstIterator dict_it = dict.Begin(); dict_it != dict_end_it; ++dict_it) { - std::string tag_key = dict_it->first; + using MetaDataStringType = itk::MetaDataObject; - // only print out if no tags have been specified, or this tag is a specified tag - if (!limit_tags || (in_tags.find(tag_key)) != in_tags.end()) + MetaDataStringType::Pointer entry = dynamic_cast(dict_it->second.GetPointer()); + + if (entry) { - std::string label_id; + std::string tag_key = dict_it->first; + + // only print out if no tags have been specified, or this tag is a specified tag + if (!limit_tags || (in_tags.find(tag_key)) != in_tags.end()) + { + std::string label_id; - const bool tag_found = itk::GDCMImageIO::GetLabelFromTag(tag_key, label_id); + const bool tag_found = itk::GDCMImageIO::GetLabelFromTag(tag_key, label_id); - std::cout << "(" << tag_key << ") " << (tag_found ? label_id.c_str() : "Unknown") << " = " << entry->GetMetaDataObjectValue() << std::endl; + std::cout << "(" << tag_key << ") " << (tag_found ? label_id.c_str() : "Unknown") << " = " << entry->GetMetaDataObjectValue() << std::endl; + } } } } + if (print_xreg_fields) + { + std::cout << "--------- xreg Basic DICOM Fields ---------" << std::endl; + + PrintDICOMFileBasicFields(ReadDICOMFileBasicFields(img_path), std::cout); + } + return kEXIT_VAL_SUCCESS; } diff --git a/docker/Dockerfile.xreg-deps b/docker/Dockerfile.xreg-deps index 7f5d6e8..673709e 100644 --- a/docker/Dockerfile.xreg-deps +++ b/docker/Dockerfile.xreg-deps @@ -18,7 +18,7 @@ RUN mkdir ffmpeg_dl && cd ffmpeg_dl && \ wget https://johnvansickle.com/ffmpeg/releases/ffmpeg-release-amd64-static.tar.xz && \ tar xf ffmpeg-release-amd64-static.tar.xz && \ mkdir -p /usr/local/bin && \ - mv ffmpeg-4.3.1-amd64-static/ffmpeg /usr/local/bin && \ + mv ffmpeg-*-amd64-static/ffmpeg /usr/local/bin && \ cd / && rm -rf ffmpeg_dl # Get a more recent version of TBB than what is probably in the package manager @@ -42,7 +42,7 @@ WORKDIR / # bring in boost 1.74 header only RUN mkdir boost_dl && cd boost_dl && \ - wget https://dl.bintray.com/boostorg/release/1.74.0/source/boost_1_74_0.tar.gz && \ + wget https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.tar.gz && \ tar xf boost_1_74_0.tar.gz && cd boost_1_74_0 && mv boost /usr/local/include && \ cd / && rm -rf boost_dl diff --git a/docker/README.md b/docker/README.md index 6db4168..d8d36a7 100644 --- a/docker/README.md +++ b/docker/README.md @@ -22,6 +22,9 @@ When these arguments are not provided, Ubuntu 16.04 is used by default. * Extracts the xReg executables for distribution. This bundles the executables and shared libraries together in an archive. ## Example Commands +A copy/pastable list of shell commands is provided in [`example_commands`](example_commands). + +Specific commands are also listed below: NOTE: all of these commands assume that the xReg repository contents are located in ~/xreg-git. These commands were tested on MacOS 10.14.6 with Docker Desktop (community) 2.3.0.5. diff --git a/docker/example_commands b/docker/example_commands new file mode 100644 index 0000000..f1be77e --- /dev/null +++ b/docker/example_commands @@ -0,0 +1,31 @@ + +# Ubuntu +export OS_NAME="ubuntu" +#export OS_VERSION="16.04" +#export OS_VERSION="18.04" +export OS_VERSION="20.04" + +# CentOS +#export OS_NAME="centos" +#export OS_VERSION="7" +#export OS_VERSION="8" + +# create the base OS image w/ compiler tools +docker build -t xreg-dev-base-${OS_NAME}-${OS_VERSION} --build-arg os_name=${OS_NAME} --build-arg os_version=${OS_VERSION} -f ~/xreg-git/docker/Dockerfile.${OS_NAME}_dev_base . + +# Grab and compile dependencies when needed +docker build -t xreg-deps-${OS_NAME}-${OS_VERSION} --build-arg os_name=${OS_NAME} --build-arg os_version=${OS_VERSION} -f ~/xreg-git/docker/Dockerfile.xreg-deps . + +# Build xReg +docker build -t xreg-${OS_NAME}-${OS_VERSION} --build-arg os_name=${OS_NAME} --build-arg os_version=${OS_VERSION} -f ~/xreg-git/docker/Dockerfile.xreg ~/xreg-git +i + +# Create xReg distributable +docker build -t xreg-dist-${OS_NAME}-${OS_VERSION} --build-arg os_name=${OS_NAME} --build-arg os_version=${OS_VERSION} -f ~/xreg-git/docker/Dockerfile.xreg-dist-bin ~/xreg-git + +# Extract the xReg distributable +XREG_DIST_CONT_ID=`docker create xreg-dist-${OS_NAME}-${OS_VERSION}` +docker cp ${XREG_DIST_CONT_ID}:xreg-${OS_NAME}-${OS_VERSION}.tar.gz . + + + diff --git a/example_build_script b/example_build_script index f6b6732..613718c 100755 --- a/example_build_script +++ b/example_build_script @@ -50,7 +50,7 @@ if [ "$NEED_TO_DOWNLOAD" = true ]; then rm $TBB_FILE - BOOST_URL="https://dl.bintray.com/boostorg/release/1.74.0/source/boost_1_74_0.zip" + BOOST_URL="https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.zip" BOOST_FILE=`basename $BOOST_URL` diff --git a/example_build_script_2 b/example_build_script_2 index 59a658f..050b7ec 100755 --- a/example_build_script_2 +++ b/example_build_script_2 @@ -88,7 +88,7 @@ if [ "$NEED_TO_DOWNLOAD" = true ]; then rm $TBB_FILE - BOOST_URL="https://dl.bintray.com/boostorg/release/1.74.0/source/boost_1_74_0.zip" + BOOST_URL="https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.zip" BOOST_FILE=`basename $BOOST_URL` diff --git a/example_build_script_win.cmd b/example_build_script_win.cmd index 828db2d..9c66321 100644 --- a/example_build_script_win.cmd +++ b/example_build_script_win.cmd @@ -28,7 +28,7 @@ REM Using a custom install version of CMake rather than the version available REM through the VS installer (3.17.20032601-MSVC_2) as the VS version resulted in REM the call to find_package() for ITK in xReg to fail. I do not know the reason, REM but using an official version downloaded from cmake.org works. -SET "PATH=C:\cmake-3.18.3-win64-x64\bin;%PATH%" +SET "PATH=C:\cmake-3.19.2-win64-x64\bin;%PATH%" REM Ninja may be installed during VS installation - it greatly speeds up the builds SET "CMAKE_GENERATOR_ARG=-G Ninja" @@ -43,6 +43,10 @@ REM I believe this is due to xReg's unconventional use of the HDF5 libraries bui REM for ITK. In the future, we should look at building HDF5 separately for use in xReg. SET "BUILD_SHARED=OFF" +REM Root directory of the xReg source code. Assumed to be located in the same +REM directory that this script was invoked from. +SET XREG_SOURCE_DIR="%~dp0" + REM create the temporary dir if necessary IF EXIST %BUILD_ROOT% ( ECHO %BUILD_ROOT% exists @@ -66,7 +70,7 @@ curl -L -O -J https://github.com/oneapi-src/oneTBB/releases/download/v2020.3/tbb tar -xf tbb-2020.3-win.zip || EXIT /b -curl -L -O -J https://dl.bintray.com/boostorg/release/1.74.0/source/boost_1_74_0.zip || EXIT /b +curl -L -O -J https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.zip || EXIT /b tar -xf boost_1_74_0.zip || EXIT /b @@ -267,12 +271,13 @@ mkdir xreg_build || EXIT /b cd xreg_build || EXIT /b -cmake %CMAKE_GENERATOR_ARG% %~dp0 ^ +cmake %CMAKE_GENERATOR_ARG% ^ -DCMAKE_PREFIX_PATH:PATH=%INSTALL_ROOT_CMAKE% ^ -DCMAKE_INSTALL_PREFIX:PATH=%INSTALL_ROOT_CMAKE% ^ -DCMAKE_CXX_STANDARD:STRING="11" ^ -DCMAKE_BUILD_TYPE:STRING=%BUILD_CONFIG% ^ - -DBUILD_SHARED_LIBS:BOOL=%BUILD_SHARED% || EXIT /b + -DBUILD_SHARED_LIBS:BOOL=%BUILD_SHARED% ^ + %XREG_SOURCE_DIR% || EXIT /b cmake --build . --config %BUILD_CONFIG% || EXIT /b diff --git a/lib/common/xregProgOptUtils.cpp b/lib/common/xregProgOptUtils.cpp index b358c8a..b1cd094 100644 --- a/lib/common/xregProgOptUtils.cpp +++ b/lib/common/xregProgOptUtils.cpp @@ -102,39 +102,35 @@ bool IsBooleanAction(const xreg::ProgOpts::StoreAction sa) } /// \brief Get the terminal dimensions in units of characters -void GetTerminalSize(xreg::ProgOpts::size_type* width, xreg::ProgOpts::size_type* height) +std::tuple +GetTerminalSize() { + // If either of the below calls fails to obtain terminal sizes, then dimensions + // of zero (e.g. no terminal) will be used. + xreg::ProgOpts::size_type width = 0; + xreg::ProgOpts::size_type height = 0; + #ifdef _WIN32 CONSOLE_SCREEN_BUFFER_INFO csbi; if (GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi)) { - *width = static_cast(csbi.dwSize.X); - *height = static_cast(csbi.dwSize.Y); - } - else - { - // TODO: determine if we have a similar case as with ioctl below, where we should - // choose some reasonable values - xregThrow("GetConsoleScreenBufferInfo() failure!"); + width = static_cast(csbi.dwSize.X); + height = static_cast(csbi.dwSize.Y); } #else struct winsize w; if (-1 != ioctl(STDOUT_FILENO, TIOCGWINSZ, &w)) { - *width = static_cast(w.ws_col); - *height = static_cast(w.ws_row); - } - else - { - // set to reasonable values - this case may occur when the -h flag is passed - // to print help, but the standard output stream is redirected to a file. - *width = 512; - *height = 512; + width = static_cast(w.ws_col); + height = static_cast(w.ws_row); } #endif + + return std::make_tuple(width, height); } -std::string FormatLines(const std::string& in, const xreg::ProgOpts::size_type indent, const xreg::ProgOpts::size_type max_width) +std::string FormatLines(const std::string& in, const xreg::ProgOpts::size_type indent, + const xreg::ProgOpts::size_type max_width) { xregASSERT(indent < max_width); // in order to print anything, we need to at least print 1 character per line @@ -198,7 +194,7 @@ void WritePrettyHelp(const xreg::ProgOpts::StringList& flag_strs, { size_type term_width = 0; size_type term_height = 0; - GetTerminalSize(&term_width, &term_height); + std::tie(term_width, term_height) = GetTerminalSize(); size_type max_flag_str_len = flag_strs[0].size(); for (size_type i = 1; i < num_flags; ++i) @@ -222,7 +218,18 @@ void WritePrettyHelp(const xreg::ProgOpts::StringList& flag_strs, } oss << desc_strs[i]; - out << FormatLines(oss.str(), max_flag_str_len, term_width); + + // Only apply formatting when the terminal width is larger than the flag length + // that we would indent by. We see terminal widths of zero on Google Colab, so it + // is important to not attempt formatting when it is not reasonable to do so. + if (max_flag_str_len < term_width) + { + out << FormatLines(oss.str(), max_flag_str_len, term_width); + } + else + { + out << oss.str() << std::endl; + } } out.flush(); diff --git a/lib/common/xregSampleUtils.cpp b/lib/common/xregSampleUtils.cpp index 7b28f0c..76de014 100644 --- a/lib/common/xregSampleUtils.cpp +++ b/lib/common/xregSampleUtils.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020, 2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -28,6 +28,8 @@ #include #include +#include "xregAssert.h" + void xreg::SeedRNGEngWithRandDev(std::mt19937* rng) { std::random_device rd; @@ -40,3 +42,141 @@ void xreg::SeedRNGEngWithRandDev(std::mt19937* rng) rng->seed(init_seed_seq); } +xreg::size_type xreg::BinCoeff(const size_type n, const size_type k) +{ + xregASSERT(k <= n); + + // https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula + + const size_type n_plus_1 = n + 1; + + const size_type limit = std::min(k, n - k); + + size_type numer = 1; + size_type denom = 1; + + for (size_type i = 1; i <= limit; ++i) + { + numer *= n_plus_1 - i; + + denom *= i; + } + + return numer / denom; +} + +std::vector> +xreg::SampleCombos(const size_type num_elem, const size_type combo_len, + const size_type num_combos, std::mt19937& rng) +{ + using Combo = std::vector; + + const size_type max_num_combos = BinCoeff(num_elem, combo_len); + + xregASSERT(num_combos <= max_num_combos); + + std::vector combos; + combos.reserve(num_combos); + + // perform a naive rejection sampling + + // this will keep track of the combos that have been sampled, and prevent them + // from being sampled again. We store a separate std::vector of the combos that + // are to be returned, as it will be free of any sorting that std::set imposes. + std::set combos_set; + + Combo tmp_combo(combo_len); + + std::uniform_int_distribution uni_dist(0, num_elem - 1); + + for (size_type combo_idx = 0; combo_idx < num_combos; ++combo_idx) + { + bool bad_combo = true; + + do + { + for (size_type i = 0; i < combo_len; ++i) + { + tmp_combo[i] = uni_dist(rng); + } + + std::sort(tmp_combo.begin(), tmp_combo.end()); + + // combo is bad if it was not inserted + bad_combo = !combos_set.insert(tmp_combo).second; + } + while (bad_combo); + + combos.push_back(tmp_combo); + } + + return combos; +} + +std::vector> +xreg::BruteForce3Combos(const size_type num_elem) +{ + xregASSERT(num_elem >= 3); + + std::vector> combos; + + std::vector tmp_arr(3); + + for (size_type i = 0; i < (num_elem - 2); ++i) + { + tmp_arr[0] = i; + + for (size_type j = (i + 1); j < (num_elem - 1); ++j) + { + tmp_arr[1] = j; + + for (size_type k = (j + 1); k < num_elem; ++k) + { + tmp_arr[2] = k; + + combos.push_back(tmp_arr); + } + } + } + + xregASSERT(combos.size() == BinCoeff(num_elem, 3)); + + return combos; +} + +std::vector> +xreg::BruteForce4Combos(const size_type num_elem) +{ + xregASSERT(num_elem >= 4); + + std::vector> combos; + + std::vector tmp_arr(4); + + for (size_type i = 0; i < (num_elem - 3); ++i) + { + tmp_arr[0] = i; + + for (size_type j = (i + 1); j < (num_elem - 2); ++j) + { + tmp_arr[1] = j; + + for (size_type k = (j + 1); k < (num_elem - 1); ++k) + { + tmp_arr[2] = k; + + for (size_type l = (k + 1); l < num_elem; ++l) + { + tmp_arr[3] = l; + + combos.push_back(tmp_arr); + } + } + } + } + + xregASSERT(combos.size() == BinCoeff(num_elem, 4)); + + return combos; +} + diff --git a/lib/common/xregSampleUtils.h b/lib/common/xregSampleUtils.h index 7653ccd..3204dc1 100644 --- a/lib/common/xregSampleUtils.h +++ b/lib/common/xregSampleUtils.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020, 2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -27,12 +27,33 @@ #include +#include "xregCommon.h" + namespace xreg { -// helper function to initialize mt rng engine with a random device, this is code that I keep having to re-write everywhere +// helper function to initialize mt rng engine with a random device, this is +// code that I keep having to re-write everywhere. void SeedRNGEngWithRandDev(std::mt19937* rng); +// Compute binomial coefficient - e.g. n choose k. +size_type BinCoeff(const size_type n, const size_type k); + +// Sample some combinations of elements +std::vector> +SampleCombos(const size_type num_elem, const size_type combo_len, + const size_type num_combos, std::mt19937& rng); + +// Return an exhaustive list of combinations of 3 elements from a collection of +// a specified list. Each combination is represented by a list of 3 indices. +std::vector> +BruteForce3Combos(const size_type num_elem); + +// Return an exhaustive list of combinations of 4 elements from a collection of +// a specified list. Each combination is represented by a list of 4 indices. +std::vector> +BruteForce4Combos(const size_type num_elem); + } // xreg #endif diff --git a/lib/file_formats/CMakeLists.txt b/lib/file_formats/CMakeLists.txt index 28d3ed4..09e0a2e 100644 --- a/lib/file_formats/CMakeLists.txt +++ b/lib/file_formats/CMakeLists.txt @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2020 Robert Grupp +# Copyright (c) 2020-2021 Robert Grupp # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -22,19 +22,29 @@ file(GLOB_RECURSE XREG_CUR_LIB_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/*.h") -add_library(xreg_file_formats OBJECT ${XREG_CUR_LIB_HEADERS} - xregDICOMUtils.cpp - xregCSVUtils.cpp - xregCIOSFusionDICOM.cpp - xregACSVUtils.cpp - xregFCSVUtils.cpp - xregSTLMeshIO.cpp - xregH5MeshIO.cpp - xregMeshIO.cpp - xregH5PAOIO.cpp - xregPAOIO.cpp - xregH5CamModelIO.cpp - xregH5ProjDataIO.cpp - xregH5SE3OptVarsIO.cpp - xregWriteVideo.cpp) +set(xreg_file_formats_src ${XREG_CUR_LIB_HEADERS} + xregDICOMUtils.cpp + xregCSVUtils.cpp + xregCIOSFusionDICOM.cpp + xregACSVUtils.cpp + xregFCSVUtils.cpp + xregSTLMeshIO.cpp + xregH5MeshIO.cpp + xregMeshIO.cpp + xregH5PAOIO.cpp + xregPAOIO.cpp + xregH5CamModelIO.cpp + xregH5ProjDataIO.cpp + xregH5SE3OptVarsIO.cpp + xregWriteVideo.cpp + xregRadRawProj.cpp + xregStaRawVol.cpp + xregReadProjDataFromDICOM.cpp) + +if (APPLE) + set(xreg_file_formats_src ${xreg_file_formats_src} + xregAppleAVFoundation.mm) +endif () + +add_library(xreg_file_formats OBJECT ${xreg_file_formats_src}) diff --git a/lib/file_formats/xregAppleAVFoundation.h b/lib/file_formats/xregAppleAVFoundation.h new file mode 100644 index 0000000..79d0472 --- /dev/null +++ b/lib/file_formats/xregAppleAVFoundation.h @@ -0,0 +1,64 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#ifndef XREGAPPLEAVFOUNDATION_H_ +#define XREGAPPLEAVFOUNDATION_H_ + +#include "xregWriteVideo.h" + +namespace xreg +{ + +class WriteImageFramesToVideoAppleAVF : public WriteImageFramesToVideo +{ +public: + void open() override; + + void close() override; + + void write(const cv::Mat& frame) override; + + ~WriteImageFramesToVideoAppleAVF() override; + +private: + void* av_asset_writer_ = nullptr; + + void* av_asset_writer_input_ = nullptr; + + void* av_assest_writer_pix_buf_adaptor_ = nullptr; + + bool input_setup_ = false; + + int num_rows_ = 0; + int num_cols_ = 0; + + int frame_type_ = 0; + + long frame_count_ = 0; +}; + +} // xreg + +#endif + diff --git a/lib/file_formats/xregAppleAVFoundation.mm b/lib/file_formats/xregAppleAVFoundation.mm new file mode 100644 index 0000000..d5133a7 --- /dev/null +++ b/lib/file_formats/xregAppleAVFoundation.mm @@ -0,0 +1,246 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregAppleAVFoundation.h" + +#import + +#include + +#include "xregAssert.h" +#include "xregExceptionUtils.h" +#include "xregFilesystemUtils.h" + +void xreg::WriteImageFramesToVideoAppleAVF::open() +{ + if (Path(dst_vid_path).exists()) + { + std::remove(dst_vid_path.c_str()); + } + + NSError* error = nil; + + AVAssetWriter* writer = [AVAssetWriter assetWriterWithURL: + [NSURL fileURLWithPath: + [NSString stringWithUTF8String:dst_vid_path.c_str()].stringByAbbreviatingWithTildeInPath + isDirectory:NO] + fileType:AVFileTypeMPEG4 + error:&error]; + + if (writer) + { + av_asset_writer_ = writer; + } + else + { + xregThrow("Failed to initialize AVAssetWriter: %s", error.localizedDescription.UTF8String); + } + + frame_count_ = 0; +} + +xreg::WriteImageFramesToVideoAppleAVF::~WriteImageFramesToVideoAppleAVF() +{ + if (input_setup_) + { + close(); + } +} + +void xreg::WriteImageFramesToVideoAppleAVF::close() +{ + AVAssetWriterInput* writer_input = static_cast(av_asset_writer_input_); + + if (writer_input) + { + [writer_input markAsFinished]; + } + else + { + xregThrow("Cannot close writer when not opened!"); + } + + AVAssetWriter* writer = static_cast(av_asset_writer_); + + if (writer) + { + __block BOOL finished_writing = NO; + + [writer finishWritingWithCompletionHandler:^{ finished_writing = YES; }]; + + while (!finished_writing) + { + [NSThread sleepForTimeInterval:0.125]; + } + + if (writer.status == AVAssetWriterStatusFailed) + { + xregThrow("Failed to finish writing with error: %s", writer.error.localizedDescription.UTF8String); + } + else if (writer.status == AVAssetWriterStatusCancelled) + { + xregThrow("Failed to finish writing - was cancelled!"); + } + else if (writer.status == AVAssetWriterStatusUnknown) + { + xregThrow("Failed to finish writing - status unknown!"); + } + else if (writer.status == AVAssetWriterStatusWriting) + { + xregThrow("Failed to finish writing - still writing!"); + } + } + else + { + xregThrow("writer object is null! cannot finish writing!"); + } + + av_asset_writer_ = nullptr; + av_asset_writer_input_ = nullptr; + av_assest_writer_pix_buf_adaptor_ = nullptr; + + input_setup_ = false; +} + +void xreg::WriteImageFramesToVideoAppleAVF::write(const cv::Mat& frame) +{ + AVAssetWriter* writer = static_cast(av_asset_writer_); + + if (writer) + { + if (!input_setup_) + { + num_rows_ = frame.rows; + num_cols_ = frame.cols; + + frame_type_ = frame.type(); + + NSNumber* frame_width = [NSNumber numberWithInt:frame.cols]; + NSNumber* frame_height = [NSNumber numberWithInt:frame.rows]; + + NSDictionary* out_settings = @{ AVVideoCodecKey:AVVideoCodecTypeH264, + AVVideoWidthKey:frame_width, + AVVideoHeightKey:frame_height }; + + + AVAssetWriterInput* writer_input = [AVAssetWriterInput + assetWriterInputWithMediaType:AVMediaTypeVideo + outputSettings:out_settings]; + xregASSERT(writer_input); + + av_asset_writer_input_ = writer_input; + + NSDictionary* src_buf_attr = @{ (__bridge NSString*) kCVPixelBufferPixelFormatTypeKey: + [NSNumber numberWithInt:((frame_type_ == CV_8UC3) ? + kCVPixelFormatType_24RGB : kCVPixelFormatType_OneComponent8)], + (__bridge NSString*) kCVPixelBufferWidthKey:frame_width, + (__bridge NSString*) kCVPixelBufferHeightKey:frame_height }; + + AVAssetWriterInputPixelBufferAdaptor* pix_buf_adaptor = [AVAssetWriterInputPixelBufferAdaptor + assetWriterInputPixelBufferAdaptorWithAssetWriterInput:writer_input + sourcePixelBufferAttributes:src_buf_attr]; + xregASSERT(pix_buf_adaptor); + + av_assest_writer_pix_buf_adaptor_ = pix_buf_adaptor; + + xregASSERT([writer canAddInput:writer_input]); + [writer addInput:writer_input]; + + if ([writer startWriting]) + { + input_setup_ = true; + + [writer startSessionAtSourceTime:kCMTimeZero]; + } + else + { + if (writer.status == AVAssetWriterStatusFailed) + { + xregThrow("Failed to start writing with error: %s", writer.error.localizedDescription.UTF8String); + } + else + { + xregThrow("Unable to start writing video (no error provided)!"); + } + } + } + + AVAssetWriterInput* writer_input = static_cast(av_asset_writer_input_); + xregASSERT(writer_input); + + AVAssetWriterInputPixelBufferAdaptor* pix_buf_adaptor = + static_cast(av_assest_writer_pix_buf_adaptor_); + + xregASSERT(pix_buf_adaptor); + xregASSERT(pix_buf_adaptor.pixelBufferPool); + + xregASSERT(num_rows_ == frame.rows); + xregASSERT(num_cols_ == frame.cols); + xregASSERT(frame_type_ == frame.type()); + + while (!writer_input.readyForMoreMediaData) + { + [NSThread sleepForTimeInterval:0.05]; + } + + CVPixelBufferRef pixel_buf = nullptr; + CVPixelBufferPoolCreatePixelBuffer(nullptr, pix_buf_adaptor.pixelBufferPool, &pixel_buf); + xregASSERT(pixel_buf); + + CVPixelBufferLockBaseAddress(pixel_buf, 0); + + cv::Mat dst_mat(num_rows_, num_cols_, frame_type_, + CVPixelBufferGetBaseAddress(pixel_buf), + CVPixelBufferGetBytesPerRow(pixel_buf)); + + if (frame.channels() == 1) + { + frame.copyTo(dst_mat); + } + else + { + xregASSERT(frame.channels() == 3); + + cv::cvtColor(frame, dst_mat, cv::COLOR_BGR2RGB); + } + + CVPixelBufferUnlockBaseAddress(pixel_buf, 0); + + // TODO: consider switching to CMTimeMakeWithSeconds + if (![pix_buf_adaptor appendPixelBuffer:pixel_buf + withPresentationTime:CMTimeMake(frame_count_, static_cast(fps))]) + { + xregThrow("Failed to append pixel buffer!"); + } + + CVPixelBufferRelease(pixel_buf); + + ++frame_count_; + } + else + { + xregThrow("cannot write a frame before AVAssetWriter setup!"); + } +} + diff --git a/lib/file_formats/xregCIOSFusionDICOM.cpp b/lib/file_formats/xregCIOSFusionDICOM.cpp index 9bf2912..cde89eb 100644 --- a/lib/file_formats/xregCIOSFusionDICOM.cpp +++ b/lib/file_formats/xregCIOSFusionDICOM.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -24,289 +24,7 @@ #include "xregCIOSFusionDICOM.h" -#include - -#include -#include - -#include "xregStringUtils.h" -#include "xregExceptionUtils.h" -#include "xregITKIOUtils.h" -#include "xregITKOpenCVUtils.h" -#include "xregOpenCVUtils.h" -#include "xregHDF5.h" - -xreg::CIOSFusionDICOMInfo xreg::ReadCIOSFusionDICOMMetadata(const std::string& path) -{ - CIOSFusionDICOMInfo meta; - - gdcm::Reader dcm_reader; - dcm_reader.SetFileName(path.c_str()); - - if (dcm_reader.CanRead()) - { - using StudyTimeAttr = gdcm::Attribute<0x0008,0x0030>; - using SeriesTimeAttr = gdcm::Attribute<0x0008,0x0031>; - using AcquisitionTimeAttr = gdcm::Attribute<0x0008,0x0032>; - - using PatNameAttr = gdcm::Attribute<0x0010,0x0010>; - using PatIdAttr = gdcm::Attribute<0x0010,0x0020>; - - using KVPAttr = gdcm::Attribute<0x0018,0x0060>; - - using DistSrcToDetAttr = gdcm::Attribute<0x0018,0x1110>; - - using TubeCurrentAttr = gdcm::Attribute<0x0018,0x1151>; - using ExposureAttr = gdcm::Attribute<0x0018,0x1152>; - using ExposureMuAsAttr = gdcm::Attribute<0x0018,0x1153>; - - using PixelSpacingAttr = gdcm::Attribute<0x0018,0x1164>; - - using FOVOriginOffAttr = gdcm::Attribute<0x0018,0x7030>; - using FOVRotAttr = gdcm::Attribute<0x0018,0x7032>; - using FOVHorizFlipAttr = gdcm::Attribute<0x0018,0x7034>; - - using GridFocalDistAttr = gdcm::Attribute<0x0018,0x704c>; - - using TubeCurrentInMuAAttr = gdcm::Attribute<0x0018,0x8151>; - - using RowsAttr = gdcm::Attribute<0x0028,0x0010>; - using ColsAttr = gdcm::Attribute<0x0028,0x0011>; - - using WinCenterAttr = gdcm::Attribute<0x0028,0x1050>; - using WinWidthAttr = gdcm::Attribute<0x0028,0x1051>; - - std::set tags_to_read; - - tags_to_read.insert(StudyTimeAttr::GetTag()); - tags_to_read.insert(SeriesTimeAttr::GetTag()); - tags_to_read.insert(AcquisitionTimeAttr::GetTag()); - - tags_to_read.insert(PatNameAttr::GetTag()); - tags_to_read.insert(PatIdAttr::GetTag()); - - tags_to_read.insert(KVPAttr::GetTag()); - - tags_to_read.insert(DistSrcToDetAttr::GetTag()); - - tags_to_read.insert(TubeCurrentAttr::GetTag()); - tags_to_read.insert(ExposureAttr::GetTag()); - tags_to_read.insert(ExposureMuAsAttr::GetTag()); - - tags_to_read.insert(PixelSpacingAttr::GetTag()); - - tags_to_read.insert(FOVOriginOffAttr::GetTag()); - tags_to_read.insert(FOVRotAttr::GetTag()); - tags_to_read.insert(FOVHorizFlipAttr::GetTag()); - - tags_to_read.insert(GridFocalDistAttr::GetTag()); - - tags_to_read.insert(TubeCurrentInMuAAttr::GetTag()); - - tags_to_read.insert(RowsAttr::GetTag()); - tags_to_read.insert(ColsAttr::GetTag()); - - tags_to_read.insert(WinCenterAttr::GetTag()); - tags_to_read.insert(WinWidthAttr::GetTag()); - - if (dcm_reader.ReadSelectedTags(tags_to_read)) - { - gdcm::DataSet& ds = dcm_reader.GetFile().GetDataSet(); - - { - StudyTimeAttr study_time_attr; - study_time_attr.SetFromDataSet(ds); - meta.study_time = StringCast(study_time_attr.GetValue()); - } - - { - SeriesTimeAttr series_time_attr; - series_time_attr.SetFromDataSet(ds); - meta.series_time = StringCast(series_time_attr.GetValue()); - } - - { - AcquisitionTimeAttr acq_time_attr; - acq_time_attr.SetFromDataSet(ds); - meta.acquisition_time = StringCast(acq_time_attr.GetValue()); - } - - { - PatNameAttr pat_name_attr; - pat_name_attr.SetFromDataSet(ds); - meta.pat_name = StringStripExtraNulls(pat_name_attr.GetValue()); - } - - { - PatIdAttr pat_id_attr; - pat_id_attr.SetFromDataSet(ds); - meta.pat_id = StringStripExtraNulls(pat_id_attr.GetValue()); - } - - { - auto de = ds.GetDataElement(KVPAttr::GetTag()); - - if (!de.IsEmpty()) - { - KVPAttr kvp_attr; - kvp_attr.SetFromDataElement(de); - - meta.kvp = kvp_attr.GetValue(); - } - } - - { - DistSrcToDetAttr dist_src_to_det_attr; - dist_src_to_det_attr.SetFromDataSet(ds); - meta.dist_src_to_det = dist_src_to_det_attr.GetValue(); - } - - { - auto de = ds.GetDataElement(TubeCurrentAttr::GetTag()); - - if (!de.IsEmpty()) - { - TubeCurrentAttr tube_current_attr; - tube_current_attr.SetFromDataElement(de); - - meta.tube_current = tube_current_attr.GetValue(); - } - } - - { - auto de = ds.GetDataElement(ExposureAttr::GetTag()); - - if (!de.IsEmpty()) - { - ExposureAttr exposure_attr; - exposure_attr.SetFromDataElement(de); - - meta.exposure = exposure_attr.GetValue(); - } - } - - { - auto de = ds.GetDataElement(ExposureMuAsAttr::GetTag()); - - if (!de.IsEmpty()) - { - ExposureMuAsAttr exposure_mu_As_attr; - exposure_mu_As_attr.SetFromDataElement(de); - - meta.exposure_mu_As = exposure_mu_As_attr.GetValue(); - } - } - - { - PixelSpacingAttr pixel_spacing_attr; - pixel_spacing_attr.SetFromDataSet(ds); - - meta.pixel_row_spacing = pixel_spacing_attr.GetValue(0); - meta.pixel_col_spacing = pixel_spacing_attr.GetValue(1); - } - - { - FOVOriginOffAttr fov_origin_off_attr; - fov_origin_off_attr.SetFromDataSet(ds); - - meta.fov_origin_row_off = fov_origin_off_attr.GetValue(0); - meta.fov_origin_col_off = fov_origin_off_attr.GetValue(1); - } - - { - FOVRotAttr fov_rot_attr; - fov_rot_attr.SetFromDataSet(ds); - - const double r = fov_rot_attr.GetValue(); - - if (r == 0) - { - meta.fov_rot = CIOSFusionDICOMInfo::kZERO; - } - else if (r == 90) - { - meta.fov_rot = CIOSFusionDICOMInfo::kNINETY; - } - else if (r == 180) - { - meta.fov_rot = CIOSFusionDICOMInfo::kONE_EIGHTY; - } - else if (r == 270) - { - meta.fov_rot = CIOSFusionDICOMInfo::kTWO_SEVENTY; - } - else - { - xregThrow("Unsupported FOV Rotation Value: %f", r); - } - } - - { - FOVHorizFlipAttr fov_horiz_flip_attr; - fov_horiz_flip_attr.SetFromDataSet(ds); - - meta.fov_horizontal_flip = fov_horiz_flip_attr.GetValue() == "YES"; - } - - { - GridFocalDistAttr grid_focal_dist_attr; - grid_focal_dist_attr.SetFromDataSet(ds); - - meta.grid_focal_dist = grid_focal_dist_attr.GetValue(); - } - - { - auto de = ds.GetDataElement(TubeCurrentInMuAAttr::GetTag()); - - if (!de.IsEmpty()) - { - TubeCurrentInMuAAttr tube_current_in_muA_attr; - tube_current_in_muA_attr.SetFromDataElement(de); - - meta.tube_current_muA = tube_current_in_muA_attr.GetValue(); - } - } - - { - RowsAttr rows_attr; - rows_attr.SetFromDataSet(ds); - - meta.rows = rows_attr.GetValue(); - } - - { - ColsAttr cols_attr; - cols_attr.SetFromDataSet(ds); - - meta.cols = cols_attr.GetValue(); - } - - { - WinCenterAttr win_center_attr; - win_center_attr.SetFromDataSet(ds); - - meta.window_center = win_center_attr.GetValue(0); - } - - { - WinWidthAttr win_width_attr; - win_width_attr.SetFromDataSet(ds); - - meta.window_width = win_width_attr.GetValue(0); - } - } - else - { - xregThrow("ReadCIOSFusionDICOMMetadata: Failed to read tags!"); - } - } - else - { - xregThrow("ReadCIOSFusionDICOMMetadata: Could not read DICOM!"); - } - - return meta; -} +#include "xregAssert.h" xreg::Mat4x4 xreg::CIOSFusionCBCTExtrins() { @@ -328,33 +46,38 @@ xreg::Mat4x4 xreg::CIOSFusionCBCTExtrins() return extrins; } -xreg::CIOSFusionDICOMInfo xreg::MakeNaiveCIOSFusionMetaDR() +xreg::DICOMFIleBasicFields xreg::MakeNaiveCIOSFusionMetaDR() { - CIOSFusionDICOMInfo meta; + DICOMFIleBasicFields meta; - meta.dist_src_to_det = 1020; + meta.dist_src_to_det_mm = 1020; - meta.pixel_row_spacing = 0.194; - meta.pixel_col_spacing = 0.194; + meta.imager_pixel_spacing = std::array{0.194, 0.194}; + + meta.row_spacing = 0.194; + meta.col_spacing = 0.194; - meta.rows = 1536; - meta.cols = 1536; + meta.num_rows = 1536; + meta.num_cols = 1536; return meta; } xreg::CameraModel -xreg::NaiveCamModelFromCIOSFusion(const CIOSFusionDICOMInfo& meta, +xreg::NaiveCamModelFromCIOSFusion(const DICOMFIleBasicFields& meta, const bool iso_center_at_origin) { + xregASSERT(bool(meta.dist_src_to_det_mm)); + xregASSERT(bool(meta.imager_pixel_spacing)); + + const auto& ps = *meta.imager_pixel_spacing; + CameraModel cam; cam.coord_frame_type = CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z; - const Mat3x3 intrins = MakeNaiveIntrins(meta.dist_src_to_det, - meta.rows, meta.cols, - meta.pixel_row_spacing, - meta.pixel_col_spacing, - true); + const Mat3x3 intrins = MakeNaiveIntrins(*meta.dist_src_to_det_mm, + meta.num_rows, meta.num_cols, + ps[1], ps[0], true); Mat4x4 extrins = Mat4x4::Identity(); @@ -363,285 +86,29 @@ xreg::NaiveCamModelFromCIOSFusion(const CIOSFusionDICOMInfo& meta, extrins = CIOSFusionCBCTExtrins(); } - cam.setup(intrins, extrins, meta.rows, meta.cols, - meta.pixel_row_spacing, meta.pixel_col_spacing); + cam.setup(intrins, extrins, meta.num_rows, meta.num_cols, ps[1], ps[0]); return cam; } xreg::CameraModel -xreg::NaiveCamModelFromCIOSFusionExtrins(const CIOSFusionDICOMInfo& meta, +xreg::NaiveCamModelFromCIOSFusionExtrins(const DICOMFIleBasicFields& meta, const Mat4x4 extrins) { + xregASSERT(bool(meta.dist_src_to_det_mm)); + xregASSERT(bool(meta.imager_pixel_spacing)); + + const auto& ps = *meta.imager_pixel_spacing; + CameraModel cam; cam.coord_frame_type = CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z; - const Mat3x3 intrins = MakeNaiveIntrins(meta.dist_src_to_det, - meta.rows, meta.cols, - meta.pixel_row_spacing, - meta.pixel_col_spacing, - true); + const Mat3x3 intrins = MakeNaiveIntrins(*meta.dist_src_to_det_mm, + meta.num_rows, meta.num_cols, + ps[1], ps[0], true); - cam.setup(intrins, extrins, meta.rows, meta.cols, - meta.pixel_row_spacing, meta.pixel_col_spacing); + cam.setup(intrins, extrins, meta.num_rows, meta.num_cols, ps[1], ps[0]); return cam; } -namespace -{ - -using namespace xreg; - -template -void ModifyImageWithCIOSFusionRotFlipFlagsHelper(const CIOSFusionDICOMInfo& meta, - itk::Image* img) -{ - using PixelType = tPixelType; - - cv::Mat img_ocv = ShallowCopyItkToOpenCV(img); - - // check for flipping, etc. - if (meta.fov_rot == CIOSFusionDICOMInfo::kNINETY) - { - xregASSERT(meta.rows == meta.cols); - - cv::Mat tmp = img_ocv.clone(); - cv::transpose(tmp, img_ocv); - FlipImageColumns(&img_ocv); - } - else if (meta.fov_rot == CIOSFusionDICOMInfo::kONE_EIGHTY) - { - FlipImageRows(&img_ocv); - FlipImageColumns(&img_ocv); - } - else if (meta.fov_rot == CIOSFusionDICOMInfo::kTWO_SEVENTY) - { - xregASSERT(meta.rows == meta.cols); - - cv::Mat tmp = img_ocv.clone(); - cv::transpose(tmp, img_ocv); - FlipImageRows(&img_ocv); - } - - // not clear if rows or columns should be flipped here - // hopefully, we do not encounter this case - if (meta.fov_horizontal_flip) - { - FlipImageColumns(&img_ocv); - } -} - -template -std::tuple::Pointer,CIOSFusionDICOMInfo> -ReadCIOSFusionDICOMHelper(const std::string& path, const bool no_rot_or_flip_based_on_meta) -{ - using PixelType = tPixelType; - using Img = itk::Image; - using ImgPtr = typename Img::Pointer; - - const CIOSFusionDICOMInfo meta = ReadCIOSFusionDICOMMetadata(path); - - ImgPtr img = ReadDICOM2DFromDisk(path); - - if (!no_rot_or_flip_based_on_meta) - { - ModifyImageWithCIOSFusionRotFlipFlagsHelper(meta, img.GetPointer()); - } - - return std::make_tuple(img, meta); -} - -template -void UpdateLandmarkMapForCIOSFusionHelper(const CIOSFusionDICOMInfo& meta, - tPointMapItr begin_it, tPointMapItr end_it, - const bool no_rot_or_flip) -{ - using PointMapIt = tPointMapItr; - - for (PointMapIt it = begin_it; it != end_it; ++it) - { - auto& p = it->second; - - // convert from physical to index - p(0) /= meta.pixel_col_spacing; - p(1) /= meta.pixel_row_spacing; - - if (!no_rot_or_flip && - (meta.fov_horizontal_flip || (meta.fov_rot != CIOSFusionDICOMInfo::kZERO))) - { - if (meta.fov_rot == CIOSFusionDICOMInfo::kNINETY) - { - std::swap(p(0), p(1)); - } - else if (meta.fov_rot == CIOSFusionDICOMInfo::kONE_EIGHTY) - { - p(0) = meta.cols - 1 - p(0); - p(1) = meta.rows - 1 - p(1); - } - else if (meta.fov_rot == CIOSFusionDICOMInfo::kTWO_SEVENTY) - { - std::swap(p(0), p(1)); - p(1) = meta.rows - 1 - p(1); - } - - // not clear if rows or columns should be flipped here - // hopefully, we do not encounter this case - if (meta.fov_horizontal_flip) - { - p(0) = meta.cols - 1 - p(0); - } - } - } -} - -} // un-named - -void xreg::ModifyImageWithCIOSFusionRotFlipFlags(const CIOSFusionDICOMInfo& meta, - itk::Image* img) -{ - ModifyImageWithCIOSFusionRotFlipFlagsHelper(meta, img); -} - -void xreg::ModifyImageWithCIOSFusionRotFlipFlags(const CIOSFusionDICOMInfo& meta, - itk::Image* img) -{ - ModifyImageWithCIOSFusionRotFlipFlagsHelper(meta, img); -} - -std::tuple::Pointer,xreg::CIOSFusionDICOMInfo> -xreg::ReadCIOSFusionDICOMUShort(const std::string& path, const bool no_rot_or_flip_based_on_meta) -{ - return ReadCIOSFusionDICOMHelper(path, no_rot_or_flip_based_on_meta); -} - -std::tuple::Pointer,xreg::CIOSFusionDICOMInfo> -xreg::ReadCIOSFusionDICOMFloat(const std::string& path, const bool no_rot_or_flip_based_on_meta) -{ - return ReadCIOSFusionDICOMHelper(path, no_rot_or_flip_based_on_meta); -} - -void xreg::UpdateLandmarkMapForCIOSFusion(const CIOSFusionDICOMInfo& meta, - LandMap2* pts, - const bool no_rot_or_flip) -{ - UpdateLandmarkMapForCIOSFusionHelper(meta, pts->begin(), pts->end(), no_rot_or_flip); -} - -void xreg::UpdateLandmarkMapForCIOSFusion(const CIOSFusionDICOMInfo& meta, - LandMap3* pts, - const bool no_rot_or_flip) -{ - UpdateLandmarkMapForCIOSFusionHelper(meta, pts->begin(), pts->end(), no_rot_or_flip); -} - - -void xreg::WriteCIOSMetaH5(const CIOSFusionDICOMInfo& meta, H5::Group* h5) -{ - WriteSingleScalarH5("study-time", meta.study_time, h5); - WriteSingleScalarH5("series-time", meta.series_time, h5); - WriteSingleScalarH5("acquisition-time", meta.acquisition_time, h5); - - WriteStringH5("pat-name", meta.pat_name, h5, false); - WriteStringH5("pat-id", meta.pat_id, h5, false); - - if (meta.kvp) - { - WriteSingleScalarH5("kvp", *meta.kvp, h5); - } - - WriteSingleScalarH5("dist-src-to-det", meta.dist_src_to_det, h5); - - if (meta.tube_current) - { - WriteSingleScalarH5("tube-current", *meta.tube_current, h5); - } - - if (meta.exposure) - { - WriteSingleScalarH5("exposure", *meta.exposure, h5); - } - - if (meta.exposure_mu_As) - { - WriteSingleScalarH5("exposure-mu-As", *meta.exposure_mu_As, h5); - } - - WriteSingleScalarH5("pixel-row-spacing", meta.pixel_row_spacing, h5); - WriteSingleScalarH5("pixel-col-spacing", meta.pixel_col_spacing, h5); - - WriteSingleScalarH5("fov-origin-row-off", meta.fov_origin_row_off, h5); - WriteSingleScalarH5("fov-origin-col-off", meta.fov_origin_col_off, h5); - - WriteSingleScalarH5("fov-rot", static_cast(meta.fov_rot), h5); - - WriteSingleScalarH5("fov-horizontal-flip", meta.fov_horizontal_flip, h5); - - WriteSingleScalarH5("grid-focal-dist", meta.grid_focal_dist, h5); - - if (meta.tube_current_muA) - { - WriteSingleScalarH5("tube-current-muA", *meta.tube_current_muA, h5); - } - - WriteSingleScalarH5("rows", meta.rows, h5); - WriteSingleScalarH5("cols", meta.cols, h5); - - WriteSingleScalarH5("window-center", meta.window_center, h5); - WriteSingleScalarH5("window-width", meta.window_width, h5); -} - -xreg::CIOSFusionDICOMInfo xreg::ReadCIOSMetaH5(const H5::Group& h5) -{ - auto read_optional_double = [&h5] (const std::string& n) - { - boost::optional to_ret; - - if (ObjectInGroupH5(n, h5)) - { - to_ret = ReadSingleScalarH5Double(n, h5); - } - - return to_ret; - }; - - CIOSFusionDICOMInfo meta; - - meta.study_time = ReadSingleScalarH5Double("study-time", h5); - meta.series_time = ReadSingleScalarH5Double("series-time", h5); - meta.acquisition_time = ReadSingleScalarH5Double("acquisition-time", h5); - - meta.pat_name = ReadStringH5("pat-name", h5); - meta.pat_id = ReadStringH5("pat-id", h5); - - meta.kvp = read_optional_double("kvp"); - - meta.dist_src_to_det = ReadSingleScalarH5Double("dist-src-to-det", h5); - - meta.tube_current = read_optional_double("tube-current"); - meta.exposure = read_optional_double("exposure"); - meta.exposure_mu_As = read_optional_double("exposure-mu-As"); - - meta.pixel_row_spacing = ReadSingleScalarH5Double("pixel-row-spacing", h5); - meta.pixel_col_spacing = ReadSingleScalarH5Double("pixel-col-spacing", h5); - - meta.fov_origin_row_off = ReadSingleScalarH5ULong("fov-origin-row-off", h5); - meta.fov_origin_col_off = ReadSingleScalarH5ULong("fov-origin-col-off", h5); - - meta.fov_rot = static_cast(ReadSingleScalarH5Int("fov-rot", h5)); - - meta.fov_horizontal_flip = ReadSingleScalarH5Bool("fov-horizontal-flip", h5); - - meta.grid_focal_dist = ReadSingleScalarH5Double("grid-focal-dist", h5); - - meta.tube_current_muA = read_optional_double("tube-current-muA"); - - meta.rows = ReadSingleScalarH5ULong("rows", h5); - meta.cols = ReadSingleScalarH5ULong("cols", h5); - - meta.window_center = ReadSingleScalarH5Double("window-center", h5); - meta.window_width = ReadSingleScalarH5Double("window-width", h5); - - return meta; -} - diff --git a/lib/file_formats/xregCIOSFusionDICOM.h b/lib/file_formats/xregCIOSFusionDICOM.h index 2fcefc3..8bbd2c1 100644 --- a/lib/file_formats/xregCIOSFusionDICOM.h +++ b/lib/file_formats/xregCIOSFusionDICOM.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -25,145 +25,36 @@ #ifndef XREGCIOSFUSIONDICOM_H_ #define XREGCIOSFUSIONDICOM_H_ -#include - -#include - -#include "xregCommon.h" +#include "xregDICOMUtils.h" #include "xregPerspectiveXform.h" -// Forward declaration -namespace H5 -{ -class Group; -} -// End Forward declaration - namespace xreg { -struct CIOSFusionDICOMInfo -{ - // these could be used for sync'ing with an external device, - // e.g. gyro, tracker, etc. - double study_time; - double series_time; - double acquisition_time; - - std::string pat_name; - std::string pat_id; - - boost::optional kvp; - - // This is the distance in mm from the source to the - // center of the detector - double dist_src_to_det; - - boost::optional tube_current; - - boost::optional exposure; - - boost::optional exposure_mu_As; - - // These are defined to be the spacings between the centers of - // pixels on the front plane of the image receptor housing. - double pixel_row_spacing; - double pixel_col_spacing; - - unsigned long fov_origin_row_off; - unsigned long fov_origin_col_off; - - enum FOVRot - { - kZERO = 0, - kNINETY = 90, - kONE_EIGHTY = 180, - kTWO_SEVENTY = 270 - }; - - FOVRot fov_rot = kZERO; - - bool fov_horizontal_flip = false; - - double grid_focal_dist; - - boost::optional tube_current_muA; - - unsigned long rows; - unsigned long cols; - - double window_center; - double window_width; -}; - -/// \brief Read CIOS Fusion DICOM metadata. -CIOSFusionDICOMInfo ReadCIOSFusionDICOMMetadata(const std::string& path); - Mat4x4 CIOSFusionCBCTExtrins(); -/// \brief Populate a CIOS Fusion metadata struct with some basic params that match -/// a typical Digital Radiograph (DR) shot. +/// \brief Populate a DICOM metadata struct with some basic params that match +/// a typical Digital Radiograph (DR) shot for the CIOS fusion. /// /// The source to detector distance, pixel spacings, and rows/columns are all set to /// nominal values. Not rotation/flipping flags are set. -CIOSFusionDICOMInfo MakeNaiveCIOSFusionMetaDR(); +DICOMFIleBasicFields MakeNaiveCIOSFusionMetaDR(); /// \brief Create a naive camera model using nominal values in the CIOS Fusion DICOM metadata. /// /// The extrinsic coordinate frame is the same as the camera coordinate frame /// (e.g. the transformation is identity), unless the second argument is passed true, in which /// case the origin will lie at the isocenter of the c-arm (previously computed offline). -CameraModel NaiveCamModelFromCIOSFusion(const CIOSFusionDICOMInfo& meta, +CameraModel NaiveCamModelFromCIOSFusion(const DICOMFIleBasicFields& meta, const bool iso_center_at_origin = false); -/// \brief Create a naive camera model using nominal values in the CIOS Fusion DICOM metadata and Extrinsic xform. +/// \brief Create a naive camera model using nominal values in the CIOS Fusion +/// DICOM metadata and Extrinsic xform. /// /// The extrinsic coordinate frame is read in from program -CameraModel NaiveCamModelFromCIOSFusionExtrins(const CIOSFusionDICOMInfo& meta, +CameraModel NaiveCamModelFromCIOSFusionExtrins(const DICOMFIleBasicFields& meta, const Mat4x4 extrins); -/// \brief Updates a 2D image with any rotating/flipping flags present -/// in a CIOS Fusion header. -/// -/// This is useful for converting label maps in the original DICOM space -/// and updating them to match a 2D/3D registration performed on the CIOS -/// Fusion DICOM. -void ModifyImageWithCIOSFusionRotFlipFlags(const CIOSFusionDICOMInfo& meta, - itk::Image* img); - -/// \brief Updates a 2D image with any rotating/flipping flags present -/// in a CIOS Fusion header. -/// -/// This is useful for converting label maps in the original DICOM space -/// and updating them to match a 2D/3D registration performed on the CIOS -/// Fusion DICOM. -void ModifyImageWithCIOSFusionRotFlipFlags(const CIOSFusionDICOMInfo& meta, - itk::Image* img); - -/// \brief Reads a CIOS Fusion DICOM from disk, returning an ITK image object and -/// another metadata structure. -std::tuple::Pointer,CIOSFusionDICOMInfo> -ReadCIOSFusionDICOMUShort(const std::string& path, const bool no_rot_or_flip_based_on_meta = false); - -/// \brief Reads a CIOS Fusion DICOM from disk, returning an ITK image object and -/// another metadata structure. -std::tuple::Pointer,CIOSFusionDICOMInfo> -ReadCIOSFusionDICOMFloat(const std::string& path, const bool no_rot_or_flip_based_on_meta = false); - -/// \brief Updates a landmark mapping from ITK physical points to continuous indices. -void UpdateLandmarkMapForCIOSFusion(const CIOSFusionDICOMInfo& meta, - LandMap2* pts, - const bool no_rot_or_flip = false); - -/// \brief Updates a landmark mapping from ITK physical points to continuous indices. -void UpdateLandmarkMapForCIOSFusion(const CIOSFusionDICOMInfo& meta, - LandMap3* pts, - const bool no_rot_or_flip = false); - -void WriteCIOSMetaH5(const CIOSFusionDICOMInfo& meta, H5::Group* h5); - -CIOSFusionDICOMInfo ReadCIOSMetaH5(const H5::Group& h5); - } // xreg #endif diff --git a/lib/file_formats/xregDICOMUtils.cpp b/lib/file_formats/xregDICOMUtils.cpp index c01b69e..5a3dd4b 100644 --- a/lib/file_formats/xregDICOMUtils.cpp +++ b/lib/file_formats/xregDICOMUtils.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -34,9 +34,11 @@ #include #include "xregAssert.h" +#include "xregHDF5.h" +#include "xregHDF5Internal.h" #include "xregStringUtils.h" -void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicFields* dcm_info) +xreg::DICOMFIleBasicFields xreg::ReadDICOMFileBasicFields(const std::string& dcm_path) { gdcm::Reader dcm_reader; dcm_reader.SetFileName(dcm_path.c_str()); @@ -45,7 +47,15 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF using PatientIDAttr = gdcm::Attribute<0x0010,0x0020>; using StudyUIDAttr = gdcm::Attribute<0x0020,0x000D>; using SeriesUIDAttr = gdcm::Attribute<0x0020,0x000E>; - using ModalityAttr = gdcm::Attribute<0x0008,0x0060>; + + using PatientNameAttr = gdcm::Attribute<0x0010,0x0010>; + + using StudyTimeAttr = gdcm::Attribute<0x0008,0x0030>; + using SeriesTimeAttr = gdcm::Attribute<0x0008,0x0031>; + using AcquisitionTimeAttr = gdcm::Attribute<0x0008,0x0032>; + using ContentTimeAttr = gdcm::Attribute<0x0008,0x0033>; + + using ModalityAttr = gdcm::Attribute<0x0008,0x0060>; using ImgPosPatAttr = gdcm::Attribute<0x0020,0x0032>; using ImgOrientPatAttr = gdcm::Attribute<0x0020,0x0037>; @@ -62,6 +72,11 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF using SeriesDescAttr = gdcm::Attribute<0x0008,0x103E>; using ImageTypeAttr = gdcm::Attribute<0x0008,0x0008>; + + using ManufacturerAttr = gdcm::Attribute<0x0008,0x0070>; + using InstitutionNameAttr = gdcm::Attribute<0x0008,0x0080>; + using InstitutionDeptAttr = gdcm::Attribute<0x0008,0x1040>; + using ManufacturerModelNameAttr = gdcm::Attribute<0x0008,0x1090>; using SecCapDevManAttr = gdcm::Attribute<0x0018,0x1016>; using SecCapDevSWVersAttr = gdcm::Attribute<0x0018,0x1019>; @@ -75,10 +90,45 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF using ConvKernelAttr = gdcm::Attribute<0x0018,0x1210>; + using BodyPartAttr = gdcm::Attribute<0x0018,0x0015>; + using ViewPosAttr = gdcm::Attribute<0x0018,0x5101>; + + using DistSrcToDetAttr = gdcm::Attribute<0x0018,0x1110>; + using DistSrcToPatAttr = gdcm::Attribute<0x0018,0x1111>; + + using KVPAttr = gdcm::Attribute<0x0018,0x0060>; + using TubeCurrentAttr = gdcm::Attribute<0x0018,0x1151>; + using ExposuremAsAttr = gdcm::Attribute<0x0018,0x1152>; + using ExposuremuAsAttr = gdcm::Attribute<0x0018,0x1153>; + using ExposureTimeAttr = gdcm::Attribute<0x0018,0x1150>; + using DoseAreaProdAttr = gdcm::Attribute<0x0018,0x115E>; + + using FOVShapeAttr = gdcm::Attribute<0x0018,0x1147>; + // We never use the FOV dims attribute, because it does not actually compile! + //using FOVDimsAttr = gdcm::Attribute<0x0018,0x1149>; + using FOVOriginOffAttr = gdcm::Attribute<0x0018,0x7030>; + using FOVRotAttr = gdcm::Attribute<0x0018,0x7032>; + using FOVHorizFlipAttr = gdcm::Attribute<0x0018,0x7034>; + + using IntensifierDiameterAttr = gdcm::Attribute<0x0018,0x1162>; + using ImagerPixelSpacingAttr = gdcm::Attribute<0x0018,0x1164>; + + using GridFocalDistAttr = gdcm::Attribute<0x0018,0x704c>; + + using WinCenterAttr = gdcm::Attribute<0x0028,0x1050>; + using WinWidthAttr = gdcm::Attribute<0x0028,0x1051>; + std::set tags_to_read; tags_to_read.insert(PatientIDAttr::GetTag()); tags_to_read.insert(StudyUIDAttr::GetTag()); tags_to_read.insert(SeriesUIDAttr::GetTag()); + tags_to_read.insert(PatientNameAttr::GetTag()); + + tags_to_read.insert(StudyTimeAttr::GetTag()); + tags_to_read.insert(SeriesTimeAttr::GetTag()); + tags_to_read.insert(AcquisitionTimeAttr::GetTag()); + tags_to_read.insert(ContentTimeAttr::GetTag()); + tags_to_read.insert(ModalityAttr::GetTag()); tags_to_read.insert(ImgPosPatAttr::GetTag()); @@ -95,6 +145,11 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF tags_to_read.insert(ImageTypeAttr::GetTag()); + tags_to_read.insert(ManufacturerAttr::GetTag()); + tags_to_read.insert(InstitutionNameAttr::GetTag()); + tags_to_read.insert(InstitutionDeptAttr::GetTag()); + tags_to_read.insert(ManufacturerModelNameAttr::GetTag()); + tags_to_read.insert(SecCapDevManAttr::GetTag()); tags_to_read.insert(SecCapDevSWVersAttr::GetTag()); @@ -106,9 +161,38 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF tags_to_read.insert(ProtoNameAttr::GetTag()); tags_to_read.insert(ConvKernelAttr::GetTag()); + + tags_to_read.insert(BodyPartAttr::GetTag()); + tags_to_read.insert(ViewPosAttr::GetTag()); + + tags_to_read.insert(DistSrcToDetAttr::GetTag()); + tags_to_read.insert(DistSrcToPatAttr::GetTag()); + + tags_to_read.insert(KVPAttr::GetTag()); + tags_to_read.insert(TubeCurrentAttr::GetTag()); + tags_to_read.insert(ExposuremAsAttr::GetTag()); + tags_to_read.insert(ExposuremuAsAttr::GetTag()); + tags_to_read.insert(ExposureTimeAttr::GetTag()); + tags_to_read.insert(DoseAreaProdAttr::GetTag()); + + tags_to_read.insert(FOVShapeAttr::GetTag()); + tags_to_read.insert(gdcm::Tag(0x0018,0x1149)); // FOV Dims + tags_to_read.insert(FOVOriginOffAttr::GetTag()); + tags_to_read.insert(FOVRotAttr::GetTag()); + tags_to_read.insert(FOVHorizFlipAttr::GetTag()); + + tags_to_read.insert(IntensifierDiameterAttr::GetTag()); + tags_to_read.insert(ImagerPixelSpacingAttr::GetTag()); + + tags_to_read.insert(GridFocalDistAttr::GetTag()); + + tags_to_read.insert(WinCenterAttr::GetTag()); + tags_to_read.insert(WinWidthAttr::GetTag()); if (dcm_reader.ReadSelectedTags(tags_to_read)) { + DICOMFIleBasicFields dcm_info; + gdcm::DataSet& ds = dcm_reader.GetFile().GetDataSet(); // using excessive scoping here just to avoid using the wrong/previous @@ -120,26 +204,71 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { PatientIDAttr pat_id_attr; pat_id_attr.SetFromDataSet(ds); - dcm_info->patient_id = StringStripExtraNulls(pat_id_attr.GetValue()); + dcm_info.patient_id = StringStripExtraNulls(pat_id_attr.GetValue()); } { StudyUIDAttr study_uid_attr; study_uid_attr.SetFromDataSet(ds); - dcm_info->study_uid = StringStripExtraNulls(study_uid_attr.GetValue()); - dcm_info->study_uid = dcm_info->study_uid.substr(0, dcm_info->study_uid.size() - 1); + dcm_info.study_uid = StringStripExtraNulls(study_uid_attr.GetValue()); + dcm_info.study_uid = dcm_info.study_uid.substr(0, dcm_info.study_uid.size() - 1); } { SeriesUIDAttr series_uid_attr; series_uid_attr.SetFromDataSet(ds); - dcm_info->series_uid = StringStripExtraNulls(series_uid_attr.GetValue()); + dcm_info.series_uid = StringStripExtraNulls(series_uid_attr.GetValue()); + } + + { + PatientNameAttr pat_name_attr; + pat_name_attr.SetFromDataSet(ds); + dcm_info.patient_name = StringStripExtraNulls(pat_name_attr.GetValue()); + } + + { + StudyTimeAttr study_time_attr; + study_time_attr.SetFromDataSet(ds); + dcm_info.study_time = StringCast(StringStripExtraNulls(study_time_attr.GetValue())); + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x0031))) + { + SeriesTimeAttr series_time_attr; + series_time_attr.SetFromDataSet(ds); + + if (series_time_attr.GetNumberOfValues() > 0) + { + dcm_info.series_time = StringCast(StringStripExtraNulls(series_time_attr.GetValue())); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x0032))) + { + AcquisitionTimeAttr aquis_time_attr; + aquis_time_attr.SetFromDataSet(ds); + + if (aquis_time_attr.GetNumberOfValues() > 0) + { + dcm_info.acquisition_time = StringCast(StringStripExtraNulls(aquis_time_attr.GetValue())); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x0033))) + { + ContentTimeAttr content_time_attr; + content_time_attr.SetFromDataSet(ds); + + if (content_time_attr.GetNumberOfValues() > 0) + { + dcm_info.content_time = StringCast(StringStripExtraNulls(content_time_attr.GetValue())); + } } { ModalityAttr modality_attr; modality_attr.SetFromDataSet(ds); - dcm_info->modality = StringStripExtraNulls(modality_attr.GetValue()); + dcm_info.modality = StringStripExtraNulls(modality_attr.GetValue()); } // get patient position, directions, and image spacing. @@ -147,60 +276,69 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF ImgPosPatAttr img_pos_pat_attr; img_pos_pat_attr.SetFromDataSet(ds); xregASSERT(img_pos_pat_attr.GetNumberOfValues() == 3); - dcm_info->img_pos_wrt_pat[0] = img_pos_pat_attr.GetValue(0); - dcm_info->img_pos_wrt_pat[1] = img_pos_pat_attr.GetValue(1); - dcm_info->img_pos_wrt_pat[2] = img_pos_pat_attr.GetValue(2); + dcm_info.img_pos_wrt_pat[0] = img_pos_pat_attr.GetValue(0); + dcm_info.img_pos_wrt_pat[1] = img_pos_pat_attr.GetValue(1); + dcm_info.img_pos_wrt_pat[2] = img_pos_pat_attr.GetValue(2); } { ImgOrientPatAttr img_orient_pat_attr; img_orient_pat_attr.SetFromDataSet(ds); xregASSERT(img_orient_pat_attr.GetNumberOfValues() == 6); - dcm_info->col_dir[0] = img_orient_pat_attr.GetValue(0); - dcm_info->col_dir[1] = img_orient_pat_attr.GetValue(1); - dcm_info->col_dir[2] = img_orient_pat_attr.GetValue(2); - dcm_info->row_dir[0] = img_orient_pat_attr.GetValue(3); - dcm_info->row_dir[1] = img_orient_pat_attr.GetValue(4); - dcm_info->row_dir[2] = img_orient_pat_attr.GetValue(5); + dcm_info.col_dir[0] = img_orient_pat_attr.GetValue(0); + dcm_info.col_dir[1] = img_orient_pat_attr.GetValue(1); + dcm_info.col_dir[2] = img_orient_pat_attr.GetValue(2); + dcm_info.row_dir[0] = img_orient_pat_attr.GetValue(3); + dcm_info.row_dir[1] = img_orient_pat_attr.GetValue(4); + dcm_info.row_dir[2] = img_orient_pat_attr.GetValue(5); } { RowsAttr rows_attr; rows_attr.SetFromDataSet(ds); xregASSERT(rows_attr.GetNumberOfValues() == 1); - dcm_info->num_rows = rows_attr.GetValue(); + dcm_info.num_rows = rows_attr.GetValue(); } { ColsAttr cols_attr; cols_attr.SetFromDataSet(ds); xregASSERT(cols_attr.GetNumberOfValues() == 1); - dcm_info->num_cols = cols_attr.GetValue(); + dcm_info.num_cols = cols_attr.GetValue(); } { PixelSpacingAttr ps_attr; ps_attr.SetFromDataSet(ds); xregASSERT(ps_attr.GetNumberOfValues() == 2); - dcm_info->col_spacing = ps_attr.GetValue(0); - dcm_info->row_spacing = ps_attr.GetValue(1); + dcm_info.col_spacing = ps_attr.GetValue(0); + dcm_info.row_spacing = ps_attr.GetValue(1); } { - dcm_info->pat_pos_valid = true; PatPosAttr pat_pos_attr; pat_pos_attr.SetFromDataSet(ds); - xregASSERT(pat_pos_attr.GetNumberOfValues() == 1); - dcm_info->pat_pos = StringStripExtraNulls(pat_pos_attr.GetValue()); + + if (pat_pos_attr.GetNumberOfValues() > 0) + { + xregASSERT(pat_pos_attr.GetNumberOfValues() == 1); + dcm_info.pat_pos = StringStripExtraNulls(pat_pos_attr.GetValue()); + } } + if (ds.FindDataElement(gdcm::Tag(0x0020,0x0020))) { - dcm_info->pat_orient_valid = true; PatOrientAttr pat_orient_attr; pat_orient_attr.SetFromDataSet(ds); - xregASSERT(pat_orient_attr.GetNumberOfValues() == 2); - dcm_info->pat_orient[0] = StringStripExtraNulls(pat_orient_attr.GetValue(0)); - dcm_info->pat_orient[1] = StringStripExtraNulls(pat_orient_attr.GetValue(1)); + + if (pat_orient_attr.GetNumberOfValues() > 0) + { + xregASSERT(pat_orient_attr.GetNumberOfValues() == 2); + + dcm_info.pat_orient = std::array { + StringStripExtraNulls(pat_orient_attr.GetValue(0)), + StringStripExtraNulls(pat_orient_attr.GetValue(1)) }; + } } { @@ -211,12 +349,7 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { xregASSERT(study_desc_attr.GetNumberOfValues() == 1); - dcm_info->study_desc_valid = true; - dcm_info->study_desc = StringStripExtraNulls(study_desc_attr.GetValue()); - } - else - { - dcm_info->study_desc_valid = false; + dcm_info.study_desc = StringStripExtraNulls(study_desc_attr.GetValue()); } } @@ -228,12 +361,7 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { xregASSERT(series_desc_attr.GetNumberOfValues() == 1); - dcm_info->series_desc_valid = true; - dcm_info->series_desc = StringStripExtraNulls(series_desc_attr.GetValue()); - } - else - { - dcm_info->series_desc_valid = false; + dcm_info.series_desc = StringStripExtraNulls(series_desc_attr.GetValue()); } } @@ -245,20 +373,69 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF if (len_img_type_attr > 0) { - dcm_info->image_type_valid = true; + std::vector image_type; - dcm_info->image_type.clear(); - dcm_info->image_type.reserve(len_img_type_attr); + image_type.reserve(len_img_type_attr); for (int img_type_idx = 0; img_type_idx < len_img_type_attr; ++img_type_idx) { - dcm_info->image_type.push_back(StringStrip(StringStripExtraNulls( + image_type.push_back(StringStrip(StringStripExtraNulls( img_type_attr.GetValue(img_type_idx)))); } + + dcm_info.image_type = image_type; + } + } + + { + ManufacturerAttr manufacturer_attr; + manufacturer_attr.SetFromDataSet(ds); + + if (manufacturer_attr.GetNumberOfValues() > 0) + { + xregASSERT(manufacturer_attr.GetNumberOfValues() == 1); + + dcm_info.manufacturer = StringStrip(StringStripExtraNulls(manufacturer_attr.GetValue())); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x0080))) + { + InstitutionNameAttr inst_name_attr; + inst_name_attr.SetFromDataSet(ds); + + if (inst_name_attr.GetNumberOfValues() > 0) + { + xregASSERT(inst_name_attr.GetNumberOfValues() == 1); + + dcm_info.institution_name = StringStrip(StringStripExtraNulls(inst_name_attr.GetValue())); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x1040))) + { + InstitutionDeptAttr inst_dept_attr; + inst_dept_attr.SetFromDataSet(ds); + + if (inst_dept_attr.GetNumberOfValues() > 0) + { + xregASSERT(inst_dept_attr.GetNumberOfValues() == 1); + + dcm_info.department_name = StringStrip(StringStripExtraNulls(inst_dept_attr.GetValue())); } - else + } + + if (ds.FindDataElement(gdcm::Tag(0x0008,0x1090))) + { + ManufacturerModelNameAttr manufacturer_model_attr; + manufacturer_model_attr.SetFromDataSet(ds); + + if (manufacturer_model_attr.GetNumberOfValues() > 0) { - dcm_info->image_type_valid = false; + xregASSERT(manufacturer_model_attr.GetNumberOfValues() == 1); + + dcm_info.manufacturers_model_name = StringStrip(StringStripExtraNulls( + manufacturer_model_attr.GetValue())); } } @@ -268,13 +445,8 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF if (sec_cap_dev_man_attr.GetNumberOfValues() > 0) { - dcm_info->sec_cap_dev_manufacturer_valid = true; xregASSERT(sec_cap_dev_man_attr.GetNumberOfValues() == 1); - dcm_info->sec_cap_dev_manufacturer = StringStripExtraNulls(sec_cap_dev_man_attr.GetValue()); - } - else - { - dcm_info->sec_cap_dev_manufacturer_valid = false; + dcm_info.sec_cap_dev_manufacturer = StringStripExtraNulls(sec_cap_dev_man_attr.GetValue()); } } @@ -284,14 +456,8 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF if (sec_cap_dev_sw_vers_attr.GetNumberOfValues() > 0) { - dcm_info->sec_cap_dev_software_versions_valid = true; - xregASSERT(sec_cap_dev_sw_vers_attr.GetNumberOfValues() == 1); - dcm_info->sec_cap_dev_software_versions = StringStripExtraNulls(sec_cap_dev_sw_vers_attr.GetValue()); - } - else - { - dcm_info->sec_cap_dev_software_versions_valid = false; + dcm_info.sec_cap_dev_software_versions = StringStripExtraNulls(sec_cap_dev_sw_vers_attr.GetValue()); } } @@ -301,25 +467,21 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF if (sw_vers_attr.GetNumberOfValues() > 0) { - dcm_info->software_versions_valid = true; - const int num_sw_ver_toks = sw_vers_attr.GetNumberOfValues(); - dcm_info->software_versions.clear(); - dcm_info->software_versions.reserve(num_sw_ver_toks); + std::vector software_versions; + software_versions.reserve(num_sw_ver_toks); for (int sw_ver_idx = 0; sw_ver_idx < num_sw_ver_toks; ++sw_ver_idx) { - dcm_info->software_versions.push_back(StringStripExtraNulls(sw_vers_attr.GetValue(sw_ver_idx))); + software_versions.push_back(StringStripExtraNulls(sw_vers_attr.GetValue(sw_ver_idx))); } - } - else - { - dcm_info->software_versions_valid = false; + + dcm_info.software_versions = software_versions; } } - // See note below about GDCM always populating fields even when they are not present + // See *** note *** below about GDCM always populating fields even when they are not present // For this case it is always set to empty, but it is more appropriate to mark as invalid if (ds.FindDataElement(gdcm::Tag(0x0008,0x9206))) { @@ -330,19 +492,11 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { xregASSERT(vol_props_attr.GetNumberOfValues() == 1); - dcm_info->vol_props_valid = true; - dcm_info->vol_props = StringStrip(StringStripExtraNulls(vol_props_attr.GetValue())); - } - else - { - dcm_info->vol_props_valid = false; + dcm_info.vol_props = StringStrip(StringStripExtraNulls(vol_props_attr.GetValue())); } } - else - { - dcm_info->vol_props_valid = false; - } - + + // *** NOTE *** // Doing things a little differently here as the call to SetFromDataSet() will always // result in a value even when the number of frames attribute is not present in the // file. In the case when it is not present, there appears to be some memory corruption @@ -356,18 +510,9 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { xregASSERT(num_frames_attr.GetNumberOfValues() == 1); - dcm_info->num_frames_valid = true; - dcm_info->num_frames = num_frames_attr.GetValue(); - } - else - { - dcm_info->num_frames_valid = false; + dcm_info.num_frames = num_frames_attr.GetValue(); } } - else - { - dcm_info->num_frames_valid = false; - } { ProtoNameAttr proto_name_attr; @@ -377,12 +522,7 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF { xregASSERT(proto_name_attr.GetNumberOfValues() == 1); - dcm_info->proto_name_valid = true; - dcm_info->proto_name = StringStripExtraNulls(proto_name_attr.GetValue()); - } - else - { - dcm_info->proto_name_valid = false; + dcm_info.proto_name = StringStripExtraNulls(proto_name_attr.GetValue()); } } @@ -392,18 +532,286 @@ void xreg::ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicF if (conv_kernel_attr.GetNumberOfValues() > 0) { - xregASSERT(conv_kernel_attr.GetNumberOfValues() == 1); + if (conv_kernel_attr.GetNumberOfValues() > 1) + { + std::cerr << "WARNING: conv. kernel attr. (0018,1210) has more than one value! " + "Only the first value will be used!" << std::endl; + } + + dcm_info.conv_kernel = StringStrip(StringStripExtraNulls(conv_kernel_attr.GetValue(0))); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x0015))) + { + BodyPartAttr body_part_attr; + body_part_attr.SetFromDataSet(ds); + + if (body_part_attr.GetNumberOfValues() > 0) + { + xregASSERT(body_part_attr.GetNumberOfValues() == 1); + + dcm_info.body_part_examined = StringStripExtraNulls(body_part_attr.GetValue()); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x5101))) + { + ViewPosAttr view_pos_attr; + view_pos_attr.SetFromDataSet(ds); + + if (view_pos_attr.GetNumberOfValues() > 0) + { + xregASSERT(view_pos_attr.GetNumberOfValues() == 1); + + dcm_info.view_position = StringStripExtraNulls(view_pos_attr.GetValue()); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1110))) + { + DistSrcToDetAttr dist_src_to_det_attr; + dist_src_to_det_attr.SetFromDataSet(ds); + + if (dist_src_to_det_attr.GetNumberOfValues() > 0) + { + xregASSERT(dist_src_to_det_attr.GetNumberOfValues() == 1); + + dcm_info.dist_src_to_det_mm = dist_src_to_det_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1111))) + { + DistSrcToPatAttr dist_src_to_pat_attr; + dist_src_to_pat_attr.SetFromDataSet(ds); + + if (dist_src_to_pat_attr.GetNumberOfValues() > 0) + { + xregASSERT(dist_src_to_pat_attr.GetNumberOfValues() == 1); + + dcm_info.dist_src_to_pat_mm = dist_src_to_pat_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x0060))) + { + KVPAttr kvp_attr; + kvp_attr.SetFromDataSet(ds); + + if (kvp_attr.GetNumberOfValues() > 0) + { + xregASSERT(kvp_attr.GetNumberOfValues() == 1); + + dcm_info.kvp = kvp_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1151))) + { + TubeCurrentAttr tube_current_attr; + tube_current_attr.SetFromDataSet(ds); + + if (tube_current_attr.GetNumberOfValues() > 0) + { + xregASSERT(tube_current_attr.GetNumberOfValues() == 1); + + dcm_info.tube_current_mA = tube_current_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1152))) + { + ExposuremAsAttr exposure_mAs_attr; + exposure_mAs_attr.SetFromDataSet(ds); + + if (exposure_mAs_attr.GetNumberOfValues() > 0) + { + xregASSERT(exposure_mAs_attr.GetNumberOfValues() == 1); + + dcm_info.exposure_mAs = exposure_mAs_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1153))) + { + ExposuremuAsAttr exposure_muAs_attr; + exposure_muAs_attr.SetFromDataSet(ds); + + if (exposure_muAs_attr.GetNumberOfValues() > 0) + { + xregASSERT(exposure_muAs_attr.GetNumberOfValues() == 1); + + dcm_info.exposure_muAs = exposure_muAs_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1150))) + { + ExposureTimeAttr exposure_time_attr; + exposure_time_attr.SetFromDataSet(ds); + + if (exposure_time_attr.GetNumberOfValues() > 0) + { + xregASSERT(exposure_time_attr.GetNumberOfValues() == 1); + + dcm_info.exposure_time_ms = exposure_time_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x115E))) + { + DoseAreaProdAttr dose_area_prod_attr; + dose_area_prod_attr.SetFromDataSet(ds); + + if (dose_area_prod_attr.GetNumberOfValues() > 0) + { + xregASSERT(dose_area_prod_attr.GetNumberOfValues() == 1); + + dcm_info.dose_area_product_dGy_cm_sq = dose_area_prod_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1147))) + { + FOVShapeAttr fov_shape_attr; + fov_shape_attr.SetFromDataSet(ds); + + if (fov_shape_attr.GetNumberOfValues() > 0) + { + xregASSERT(fov_shape_attr.GetNumberOfValues() == 1); + + dcm_info.fov_shape = StringStrip(StringStripExtraNulls(fov_shape_attr.GetValue())); + } + } + + { + // See note above why we do things differently for FOV dims + gdcm::Tag fov_dims_tag(0x0018,0x1149); + + if (ds.FindDataElement(fov_dims_tag)) + { + // The FOV dims value is stored as an integer string or two integer strings, + // separated by \, when the FOV shape is RECTANGLE + + try + { + std::vector val_str = dynamic_cast( + ds.GetDataElement(fov_dims_tag).GetValue()); + val_str.push_back(0); + + dcm_info.fov_dims = StringCast(StringSplit(val_str.data(), "\\")); + } + catch (const std::bad_cast&) { } + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x7030))) + { + FOVOriginOffAttr fov_origin_off_attr; + fov_origin_off_attr.SetFromDataSet(ds); + + if (fov_origin_off_attr.GetNumberOfValues() > 0) + { + xregASSERT(fov_origin_off_attr.GetNumberOfValues() == 2); + + std::array tmp_fov_origin_off; + + tmp_fov_origin_off[0] = fov_origin_off_attr.GetValue(0); + tmp_fov_origin_off[1] = fov_origin_off_attr.GetValue(1); + + dcm_info.fov_origin_off = tmp_fov_origin_off; + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x7032))) + { + FOVRotAttr fov_rot_attr; + fov_rot_attr.SetFromDataSet(ds); + + if (fov_rot_attr.GetNumberOfValues() > 0) + { + const int tmp_fov_rot = fov_rot_attr.GetValue(); + + xregASSERT((tmp_fov_rot == 0) || (tmp_fov_rot == 90) || + (tmp_fov_rot == 180) || (tmp_fov_rot == 270)); + + dcm_info.fov_rot = static_cast(tmp_fov_rot); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x7034))) + { + FOVHorizFlipAttr fov_horiz_flip_attr; + fov_horiz_flip_attr.SetFromDataSet(ds); + + if (fov_horiz_flip_attr.GetNumberOfValues() > 0) + { + dcm_info.fov_horizontal_flip = fov_horiz_flip_attr.GetValue() == "YES"; + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1162))) + { + IntensifierDiameterAttr intensifier_diam_attr; + intensifier_diam_attr.SetFromDataSet(ds); + + if (intensifier_diam_attr.GetNumberOfValues() > 0) + { + xregASSERT(intensifier_diam_attr.GetNumberOfValues() == 1); + + dcm_info.intensifier_diameter_mm = intensifier_diam_attr.GetValue(); + } + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x1164))) + { + ImagerPixelSpacingAttr imager_pixel_spacing_attr; + imager_pixel_spacing_attr.SetFromDataSet(ds); + + xregASSERT(imager_pixel_spacing_attr.GetNumberOfValues() == 2); + + dcm_info.imager_pixel_spacing = + std::array { static_cast(imager_pixel_spacing_attr.GetValue(0)), + static_cast(imager_pixel_spacing_attr.GetValue(1)) }; + } + + if (ds.FindDataElement(gdcm::Tag(0x0018,0x704c))) + { + GridFocalDistAttr grid_focal_dist_attr; + grid_focal_dist_attr.SetFromDataSet(ds); - dcm_info->conv_kernel_valid = true; - dcm_info->conv_kernel = StringStrip(StringStripExtraNulls(conv_kernel_attr.GetValue())); + if (grid_focal_dist_attr.GetNumberOfValues() > 0) + { + dcm_info.grid_focal_dist_mm = grid_focal_dist_attr.GetValue(); } - else + } + + if (ds.FindDataElement(gdcm::Tag(0x0028,0x1050))) + { + WinCenterAttr win_center_attr; + win_center_attr.SetFromDataSet(ds); + + if (win_center_attr.GetNumberOfValues() > 0) { - dcm_info->conv_kernel_valid = false; + dcm_info.window_center = win_center_attr.GetValue(0); } } - dcm_info->file_path = dcm_path; + if (ds.FindDataElement(gdcm::Tag(0x0028,0x1051))) + { + WinWidthAttr win_width_attr; + win_width_attr.SetFromDataSet(ds); + + if (win_width_attr.GetNumberOfValues() > 0) + { + dcm_info.window_width = win_width_attr.GetValue(0); + } + } + + dcm_info.file_path = dcm_path; + + return dcm_info; } else { @@ -421,57 +829,111 @@ void xreg::PrintDICOMFileBasicFields(const DICOMFIleBasicFields& dcm_info, std:: { const std::string kNOT_PROVIDED_STR(""); - out << indent << " File Path: " << dcm_info.file_path << '\n' - << indent << " Patient ID: " << dcm_info.patient_id << '\n' - << indent << " Study UID: " << dcm_info.study_uid << '\n' - << indent << " Series UID: " << dcm_info.series_uid << '\n' - << indent << " Modality: " << dcm_info.modality << '\n' - << indent << " Image Position: " << fmt::sprintf("[%+10.4f,%+10.4f,%+10.4f]", dcm_info.img_pos_wrt_pat[0], dcm_info.img_pos_wrt_pat[1], dcm_info.img_pos_wrt_pat[2]) << '\n' - << indent << " Image Col Dir.: " << fmt::sprintf("[%+0.4f,%+0.4f,%+0.4f]", dcm_info.col_dir[0], dcm_info.col_dir[1], dcm_info.col_dir[2]) << '\n' - << indent << " Image Row Dir.: " << fmt::sprintf("[%+0.4f,%+0.4f,%+0.4f]", dcm_info.row_dir[0], dcm_info.row_dir[1], dcm_info.row_dir[2]) << '\n' - << indent << " Image Col Spacing: " << fmt::sprintf("%0.4f", dcm_info.col_spacing) << '\n' - << indent << " Image Row Spacing: " << fmt::sprintf("%0.4f", dcm_info.row_spacing) << '\n' - << indent << " Image Num Rows: " << dcm_info.num_rows << '\n' - << indent << " Image Num Cols: " << dcm_info.num_cols << '\n' - << indent << " Patient Position: " << (dcm_info.pat_pos_valid ? dcm_info.pat_pos : kNOT_PROVIDED_STR) << '\n' - << indent << " Patient Orient.: " << (dcm_info.pat_orient_valid ? fmt::sprintf("%s , %s", dcm_info.pat_orient[0], dcm_info.pat_orient[1]) : kNOT_PROVIDED_STR) << '\n' - << indent << " Study Desc.: " << (dcm_info.study_desc_valid ? dcm_info.study_desc : kNOT_PROVIDED_STR) << '\n' - << indent << " Series Desc.: " << (dcm_info.series_desc_valid ? dcm_info.series_desc : kNOT_PROVIDED_STR) << '\n' - << indent << " Image Type: " << (dcm_info.image_type_valid ? JoinTokens(dcm_info.image_type, " , ") : kNOT_PROVIDED_STR) << '\n' - << indent << " Sec. Cap. Dev. Man.: " << (dcm_info.sec_cap_dev_manufacturer_valid ? dcm_info.sec_cap_dev_manufacturer : kNOT_PROVIDED_STR) << '\n' - << indent << " Sec. Cap. Dev. SW Ver.: " << (dcm_info.sec_cap_dev_software_versions_valid ? dcm_info.sec_cap_dev_software_versions : kNOT_PROVIDED_STR) << '\n' - << indent << " Software Versions: " << (dcm_info.software_versions_valid ? JoinTokens(dcm_info.software_versions, " , ") : kNOT_PROVIDED_STR) << '\n' - << indent << " Vol. Props.: " << (dcm_info.vol_props_valid ? dcm_info.vol_props : kNOT_PROVIDED_STR) << '\n' - << indent << " Num. Frames: " << (dcm_info.num_frames_valid ? fmt::format("{}", dcm_info.num_frames) : kNOT_PROVIDED_STR) << '\n' - << indent << " Protocol Name: " << (dcm_info.proto_name_valid ? dcm_info.proto_name : kNOT_PROVIDED_STR) << '\n' - << indent << " Conv. Kernel: " << (dcm_info.conv_kernel_valid ? dcm_info.conv_kernel : kNOT_PROVIDED_STR) << '\n'; + out << indent << " File Path: " << dcm_info.file_path << '\n' + << indent << " Patient ID: " << dcm_info.patient_id << '\n' + << indent << " Patient Name: " << dcm_info.patient_name << '\n' + << indent << " Study UID: " << dcm_info.study_uid << '\n' + << indent << " Series UID: " << dcm_info.series_uid << '\n' + << indent << " Study Time: " << dcm_info.study_time << '\n' + << indent << " Series Time: " << (dcm_info.series_time ? fmt::format("{}", *dcm_info.series_time) : kNOT_PROVIDED_STR) << '\n' + << indent << " Acquisition Time: " << (dcm_info.acquisition_time ? fmt::format("{}", *dcm_info.acquisition_time) : kNOT_PROVIDED_STR) << '\n' + << indent << " Content Time: " << (dcm_info.content_time ? fmt::format("{}", *dcm_info.content_time) : kNOT_PROVIDED_STR) << '\n' + << indent << " Modality: " << dcm_info.modality << '\n' + << indent << " Image Position: " << fmt::sprintf("[%+10.4f,%+10.4f,%+10.4f]", dcm_info.img_pos_wrt_pat[0], dcm_info.img_pos_wrt_pat[1], dcm_info.img_pos_wrt_pat[2]) << '\n' + << indent << " Image Col Dir.: " << fmt::sprintf("[%+0.4f,%+0.4f,%+0.4f]", dcm_info.col_dir[0], dcm_info.col_dir[1], dcm_info.col_dir[2]) << '\n' + << indent << " Image Row Dir.: " << fmt::sprintf("[%+0.4f,%+0.4f,%+0.4f]", dcm_info.row_dir[0], dcm_info.row_dir[1], dcm_info.row_dir[2]) << '\n' + << indent << " Image Col Spacing: " << fmt::sprintf("%0.4f", dcm_info.col_spacing) << '\n' + << indent << " Image Row Spacing: " << fmt::sprintf("%0.4f", dcm_info.row_spacing) << '\n' + << indent << " Image Num Rows: " << dcm_info.num_rows << '\n' + << indent << " Image Num Cols: " << dcm_info.num_cols << '\n' + << indent << " Patient Position: " << (dcm_info.pat_pos ? *dcm_info.pat_pos : kNOT_PROVIDED_STR) << '\n' + << indent << " Patient Orient.: " << (dcm_info.pat_orient ? + fmt::sprintf("%s , %s", + (*dcm_info.pat_orient)[0], + (*dcm_info.pat_orient)[1]) : + kNOT_PROVIDED_STR) << '\n' + << indent << " Study Desc.: " << (dcm_info.study_desc ? *dcm_info.study_desc : kNOT_PROVIDED_STR) << '\n' + << indent << " Series Desc.: " << (dcm_info.series_desc ? *dcm_info.series_desc : kNOT_PROVIDED_STR) << '\n' + << indent << " Image Type: " << (dcm_info.image_type ? JoinTokens(*dcm_info.image_type, " , ") : kNOT_PROVIDED_STR) << '\n' + << indent << " Manufacturer: " << dcm_info.manufacturer << '\n' + << indent << " Institution Name: " << (dcm_info.institution_name ? *dcm_info.institution_name : kNOT_PROVIDED_STR) << '\n' + << indent << " Department Name: " << (dcm_info.department_name ? *dcm_info.department_name : kNOT_PROVIDED_STR) << '\n' + << indent << " Manuf. Model Name: " << (dcm_info.manufacturers_model_name ? *dcm_info.manufacturers_model_name : kNOT_PROVIDED_STR) << '\n' + << indent << " Sec. Cap. Dev. Man.: " << (dcm_info.sec_cap_dev_manufacturer ? *dcm_info.sec_cap_dev_manufacturer : kNOT_PROVIDED_STR) << '\n' + << indent << " Sec. Cap. Dev. SW Ver.: " << (dcm_info.sec_cap_dev_software_versions ? *dcm_info.sec_cap_dev_software_versions : kNOT_PROVIDED_STR) << '\n' + << indent << " Software Versions: " << (dcm_info.software_versions ? JoinTokens(*dcm_info.software_versions, " , ") : kNOT_PROVIDED_STR) << '\n' + << indent << " Vol. Props.: " << (dcm_info.vol_props ? *dcm_info.vol_props : kNOT_PROVIDED_STR) << '\n' + << indent << " Num. Frames: " << (dcm_info.num_frames ? fmt::format("{}", *dcm_info.num_frames) : kNOT_PROVIDED_STR) << '\n' + << indent << " Protocol Name: " << (dcm_info.proto_name ? *dcm_info.proto_name : kNOT_PROVIDED_STR) << '\n' + << indent << " Conv. Kernel: " << (dcm_info.conv_kernel ? *dcm_info.conv_kernel : kNOT_PROVIDED_STR) << '\n' + << indent << " Body Part: " << (dcm_info.body_part_examined ? *dcm_info.body_part_examined : kNOT_PROVIDED_STR) << '\n' + << indent << " View Position: " << (dcm_info.view_position ? *dcm_info.view_position : kNOT_PROVIDED_STR) << '\n' + << indent << " Dist. Src-to-Det. (mm): " << (dcm_info.dist_src_to_det_mm ? fmt::format("{:.1f}", *dcm_info.dist_src_to_det_mm) : kNOT_PROVIDED_STR) << '\n' + << indent << " Dist. Src-to-Pat. (mm): " << (dcm_info.dist_src_to_pat_mm ? fmt::format("{:.1f}", *dcm_info.dist_src_to_pat_mm) : kNOT_PROVIDED_STR) << '\n' + << indent << " Peak Engergy (kVp): " << (dcm_info.kvp ? fmt::format("{:.1f}", *dcm_info.kvp) : kNOT_PROVIDED_STR) << '\n' + << indent << " Tube Current (mA): " << (dcm_info.tube_current_mA ? fmt::format("{:.1f}", *dcm_info.tube_current_mA) : kNOT_PROVIDED_STR) << '\n' + << indent << " Exposure (mAs): " << (dcm_info.exposure_mAs ? fmt::format("{:.3f}", *dcm_info.exposure_mAs) : kNOT_PROVIDED_STR) << '\n' + << indent << " Exposure (muAs): " << (dcm_info.exposure_muAs ? fmt::format("{:.3f}", *dcm_info.exposure_muAs) : kNOT_PROVIDED_STR) << '\n' + << indent << " Exposure Time (ms): " << (dcm_info.exposure_time_ms ? fmt::format("{:.3f}", *dcm_info.exposure_time_ms) : kNOT_PROVIDED_STR) << '\n' + << indent << " Dose Area Prod. (dGy*cm^2): " << (dcm_info.dose_area_product_dGy_cm_sq ? fmt::format("{:.3f}", *dcm_info.dose_area_product_dGy_cm_sq) : kNOT_PROVIDED_STR) << '\n' + << indent << " FOV Shape: " << (dcm_info.fov_shape ? *dcm_info.fov_shape : kNOT_PROVIDED_STR) << '\n' + << indent << " FOV Dimensions (mm): " << (dcm_info.fov_dims ? JoinTokens(ToStrings(*dcm_info.fov_dims), " , ") : kNOT_PROVIDED_STR) << '\n' + << indent << " FOV Origin Offset: " << (dcm_info.fov_origin_off ? fmt::format("[{} , {}]", (*dcm_info.fov_origin_off)[0], (*dcm_info.fov_origin_off)[1]) : kNOT_PROVIDED_STR) << '\n' + << indent << " FOV Rotation: " << (dcm_info.fov_rot ? fmt::format("{}", static_cast(*dcm_info.fov_rot)) : kNOT_PROVIDED_STR) << '\n' + << indent << " FOV Horizontal Flip: " << (dcm_info.fov_horizontal_flip ? std::string(*dcm_info.fov_horizontal_flip ? "YES" : "NO") : kNOT_PROVIDED_STR) << '\n' + << indent << " Intensifier Diam. (mm): " << (dcm_info.intensifier_diameter_mm ? fmt::format("{:.1f}", *dcm_info.intensifier_diameter_mm) : kNOT_PROVIDED_STR) << '\n' + << indent << " Imager Pix. Spacing (mm/px): " + << (dcm_info.imager_pixel_spacing ? + fmt::format("[ {:.3f} , {:.3f}]", + (*dcm_info.imager_pixel_spacing)[0], + (*dcm_info.imager_pixel_spacing)[1]) : + kNOT_PROVIDED_STR) + << '\n' + << indent << " Grid Focal Dist. (mm): " << (dcm_info.grid_focal_dist_mm ? fmt::format("{}", *dcm_info.grid_focal_dist_mm) : kNOT_PROVIDED_STR) << '\n' + << indent << " Window Center: " << (dcm_info.window_center ? fmt::format("{}", *dcm_info.window_center) : kNOT_PROVIDED_STR) << '\n' + << indent << " Window Width: " << (dcm_info.window_width ? fmt::format("{}", *dcm_info.window_width) : kNOT_PROVIDED_STR) << '\n' + << '\n'; out.flush(); } bool xreg::IsLocalizer(const DICOMFIleBasicFields& dcm_info) { - return dcm_info.image_type_valid && - (std::find(dcm_info.image_type.begin(), dcm_info.image_type.end(), "LOCALIZER") - != dcm_info.image_type.end()); + return dcm_info.image_type && + (std::find(dcm_info.image_type->begin(), dcm_info.image_type->end(), "LOCALIZER") + != dcm_info.image_type->end()); } bool xreg::IsMRLocalizer(const DICOMFIleBasicFields& dcm_info) { return (dcm_info.modality == "MR") && - (dcm_info.image_type_valid && - (std::find(dcm_info.image_type.begin(), dcm_info.image_type.end(), "LOCALIZER") - != dcm_info.image_type.end())); + (dcm_info.image_type && + (std::find(dcm_info.image_type->begin(), dcm_info.image_type->end(), "LOCALIZER") + != dcm_info.image_type->end())); } bool xreg::IsVolDICOMFile(const DICOMFIleBasicFields& dcm_info) { - return dcm_info.vol_props_valid && (dcm_info.vol_props == "VOLUME"); + return dcm_info.vol_props && (*dcm_info.vol_props == "VOLUME"); } bool xreg::IsMultiFrameDICOMFile(const DICOMFIleBasicFields& dcm_info) { - return dcm_info.num_frames_valid && (dcm_info.num_frames > 1); + return dcm_info.num_frames && (*dcm_info.num_frames > 1); +} + +bool xreg::IsSecondaryDICOMFile(const DICOMFIleBasicFields& dcm_info) +{ + return dcm_info.image_type && + (std::find(dcm_info.image_type->begin(), dcm_info.image_type->end(), "SECONDARY") + != dcm_info.image_type->end()); +} + +bool xreg::IsDerivedDICOMFile(const DICOMFIleBasicFields& dcm_info) +{ + return dcm_info.image_type && + (std::find(dcm_info.image_type->begin(), dcm_info.image_type->end(), "DERIVED") + != dcm_info.image_type->end()); } void xreg::GetDICOMDirs(const std::string& root_dir_path, PathStringList* dir_paths) @@ -572,8 +1034,13 @@ void xreg::GetDICOMFilePathObjsInDir(const std::string& dir, PathList* dcm_paths void xreg::GetOrgainizedDICOMInfos(const std::string& root_dir_path, OrganizedDICOMFiles* org_dcm, const bool inc_localizer, - const bool inc_multi_frame_files) + const bool inc_multi_frame_files, + const bool inc_secondary, + const bool inc_derived, + const std::vector& modalities) { + const bool check_modality = !modalities.empty(); + org_dcm->patient_infos.clear(); org_dcm->root_dir = root_dir_path; @@ -593,10 +1060,15 @@ void xreg::GetOrgainizedDICOMInfos(const std::string& root_dir_path, { const std::string cur_file_path = dcm_file_path.string(); - ReadDICOMFileBasicFields(cur_file_path, &tmp_basic_fields); + tmp_basic_fields = ReadDICOMFileBasicFields(cur_file_path); if ((inc_localizer || !IsLocalizer(tmp_basic_fields)) && - (inc_multi_frame_files || !IsMultiFrameDICOMFile(tmp_basic_fields))) + (inc_multi_frame_files || !IsMultiFrameDICOMFile(tmp_basic_fields)) && + (inc_secondary || !IsSecondaryDICOMFile(tmp_basic_fields)) && + (inc_derived || !IsDerivedDICOMFile(tmp_basic_fields)) && + (!check_modality || + (std::find(modalities.begin(), modalities.end(), tmp_basic_fields.modality) + != modalities.end()))) { org_dcm->patient_infos[tmp_basic_fields.patient_id] [tmp_basic_fields.study_uid] @@ -644,13 +1116,13 @@ void xreg::PrintOrganizedDICOMFiles(const OrganizedDICOMFiles& org_dcm, std::ost if (print_dcm_info) { - ReadDICOMFileBasicFields(paths[0], &dcm_info); + dcm_info = ReadDICOMFileBasicFields(paths[0]); if (first_series) { out << indent3 << "Study Description: " - << (dcm_info.study_desc_valid ? dcm_info.study_desc.c_str() - : "(Not Provided)") + << (dcm_info.study_desc ? dcm_info.study_desc->c_str() + : "(Not Provided)") << std::endl; } } @@ -661,8 +1133,8 @@ void xreg::PrintOrganizedDICOMFiles(const OrganizedDICOMFiles& org_dcm, std::ost if (print_dcm_info) { out << indent4 << "Series Description: " - << (dcm_info.series_desc_valid ? dcm_info.series_desc.c_str() - : "(Not Provided)") + << (dcm_info.series_desc ? dcm_info.series_desc->c_str() + : "(Not Provided)") << std::endl; } @@ -694,12 +1166,12 @@ void xreg::ReadDICOMInfosFromDir(const std::string& dir_path, DICOMFIleBasicFiel for (size_type dcm_idx = 0; dcm_idx < num_dcm; ++dcm_idx) { - ReadDICOMFileBasicFields(dcm_paths[dcm_idx], &dcms[dcm_idx]); + dcms[dcm_idx] = ReadDICOMFileBasicFields(dcm_paths[dcm_idx]); } } bool xreg::ReorderAndCheckDICOMInfos::operator()(const DICOMFIleBasicFieldsList& src_infos, - DICOMFIleBasicFieldsList* dst_infos) + DICOMFIleBasicFieldsList* dst_infos) { const CoordScalar kTOL = 1.0e-6; @@ -892,3 +1364,358 @@ bool xreg::ReorderAndCheckDICOMInfos::operator()(const DICOMFIleBasicFieldsList& return single_out_of_plane_axis && const_in_plane_spacings && const_in_plane_dims; } + +namespace +{ + +template +void WriteOptionalScalarH5(const std::string& field_name, + const boost::optional& opt_val, + H5::Group* h5) +{ + if (opt_val) + { + xreg::WriteSingleScalarH5(field_name, *opt_val, h5); + } +} + +void WriteOptionalStringH5(const std::string& field_name, + const boost::optional& opt_str, + H5::Group* h5) +{ + if (opt_str) + { + xreg::WriteStringH5(field_name, *opt_str, h5); + } +} + +void WriteOptionalListOfStrings(const std::string& field_name, + const boost::optional>& opt_strs, + H5::Group* h5) +{ + using namespace xreg; + + if (opt_strs) + { + H5::Group strs_g = h5->createGroup(field_name); + + const size_type len = opt_strs->size(); + + SetScalarAttr("len", static_cast(len), &strs_g); + + for (size_type i = 0; i < len; ++i) + { + WriteStringH5(fmt::format("{:02d}", i), opt_strs->operator[](i), &strs_g); + } + } +} + +template +boost::optional ReadOptionalScalarH5(const std::string& field_name, + const H5::Group& h5) +{ + using Scalar = tScalar; + + boost::optional opt_val; + + if (xreg::ObjectInGroupH5(field_name, h5)) + { + opt_val = xreg::detail::ReadSingleScalarH5Helper(field_name, h5); + } + + return opt_val; +} + +boost::optional ReadOptionalStringH5(const std::string& field_name, + const H5::Group& h5) +{ + boost::optional opt_str; + + if (xreg::ObjectInGroupH5(field_name, h5)) + { + opt_str = xreg::ReadStringH5(field_name, h5); + } + + return opt_str; +} + +boost::optional> +ReadOptionalListOfStrings(const std::string& field_name, + const H5::Group& h5) +{ + using namespace xreg; + + boost::optional> opt_strs; + + if (ObjectInGroupH5(field_name, h5)) + { + H5::Group strs_g = h5.openGroup(field_name); + + const size_type len = GetScalarLongAttr("len", strs_g); + + std::vector strs; + strs.reserve(len); + + for (size_type i = 0; i < len; ++i) + { + strs.push_back(ReadStringH5(fmt::format("{:02d}", i), strs_g)); + } + + opt_strs = strs; + } + + return opt_strs; +} + +} // un-named + +void xreg::WriteDICOMFieldsH5(const DICOMFIleBasicFields& dcm_info, H5::Group* h5) +{ + WriteStringH5("file-path", dcm_info.file_path, h5); + + WriteStringH5("patient-id", dcm_info.patient_id, h5); + WriteStringH5("series-uid", dcm_info.series_uid, h5); + WriteStringH5("study-uid", dcm_info.study_uid, h5); + + WriteStringH5("patient-name", dcm_info.patient_name, h5); + + WriteSingleScalarH5("study-time", dcm_info.study_time, h5); + WriteOptionalScalarH5("series-time", dcm_info.series_time, h5); + WriteOptionalScalarH5("acquisition-time", dcm_info.acquisition_time, h5); + WriteOptionalScalarH5("content-time", dcm_info.content_time, h5); + + WriteStringH5("modality", dcm_info.modality, h5); + + WriteMatrixH5("img-pos-wrt-pat", dcm_info.img_pos_wrt_pat, h5); + + WriteMatrixH5("row-dir", dcm_info.row_dir, h5); + WriteMatrixH5("col-dir", dcm_info.col_dir, h5); + + WriteSingleScalarH5("row-spacing", dcm_info.row_spacing, h5); + WriteSingleScalarH5("col-spacing", dcm_info.col_spacing, h5); + + WriteSingleScalarH5("num-rows", dcm_info.num_rows, h5); + WriteSingleScalarH5("num-cols", dcm_info.num_cols, h5); + + WriteOptionalStringH5("pat-pos", dcm_info.pat_pos, h5); + + if (dcm_info.pat_orient) + { + const auto& pat_orient = *dcm_info.pat_orient; + + WriteStringH5("pat-orient-x", pat_orient[0], h5); + WriteStringH5("pat-orient-y", pat_orient[1], h5); + } + + WriteOptionalStringH5("study-desc", dcm_info.study_desc, h5); + + WriteOptionalStringH5("series-desc", dcm_info.series_desc, h5); + + WriteOptionalListOfStrings("image-type", dcm_info.image_type, h5); + + WriteStringH5("manufacturer", dcm_info.manufacturer, h5); + + WriteOptionalStringH5("institution-name", dcm_info.institution_name, h5); + + WriteOptionalStringH5("department-name", dcm_info.department_name, h5); + + WriteOptionalStringH5("manufacturers-model-name", dcm_info.manufacturers_model_name, h5); + + WriteOptionalStringH5("sec-cap-dev-manufacturer", dcm_info.sec_cap_dev_manufacturer, h5); + + WriteOptionalStringH5("sec-cap-dev-software-versions", dcm_info.sec_cap_dev_software_versions, h5); + + WriteOptionalListOfStrings("software-versions", dcm_info.software_versions, h5); + + WriteOptionalStringH5("vol-props", dcm_info.vol_props, h5); + + WriteOptionalScalarH5("num-frames", dcm_info.num_frames, h5); + + WriteOptionalStringH5("proto-name", dcm_info.proto_name, h5); + + WriteOptionalStringH5("conv-kernel", dcm_info.conv_kernel, h5); + + WriteOptionalStringH5("body-part-examined", dcm_info.body_part_examined, h5); + + WriteOptionalStringH5("view-position", dcm_info.view_position, h5); + + WriteOptionalScalarH5("dist-src-to-det-mm", dcm_info.dist_src_to_det_mm, h5); + + WriteOptionalScalarH5("dist-src-to-pat-mm", dcm_info.dist_src_to_pat_mm, h5); + + WriteOptionalScalarH5("kvp", dcm_info.kvp, h5); + + WriteOptionalScalarH5("tube-current-mA", dcm_info.tube_current_mA, h5); + + WriteOptionalScalarH5("exposure-mAs", dcm_info.exposure_mAs, h5); + + WriteOptionalScalarH5("exposure-muAs", dcm_info.exposure_muAs, h5); + + WriteOptionalScalarH5("exposure-time-ms", dcm_info.exposure_time_ms, h5); + + WriteOptionalScalarH5("dose-area-product-dGy-cm-sq", dcm_info.dose_area_product_dGy_cm_sq, h5); + + WriteOptionalStringH5("fov-shape", dcm_info.fov_shape, h5); + + if (dcm_info.fov_dims) + { + WriteVectorH5("fov-dims", *dcm_info.fov_dims, h5); + } + + if (dcm_info.fov_origin_off) + { + const auto origin_off = *dcm_info.fov_origin_off; + + WriteSingleScalarH5("fov-origin-off-rows", origin_off[0], h5); + WriteSingleScalarH5("fov-origin-off-cols", origin_off[1], h5); + } + + if (dcm_info.fov_rot) + { + WriteSingleScalarH5("fov-rot", static_cast(*dcm_info.fov_rot), h5); + } + + WriteOptionalScalarH5("fov-horizontal-flip", dcm_info.fov_horizontal_flip, h5); + + WriteOptionalScalarH5("intensifier-diameter-mm", dcm_info.intensifier_diameter_mm, h5); + + if (dcm_info.imager_pixel_spacing) + { + const auto& ps = *dcm_info.imager_pixel_spacing; + + WriteSingleScalarH5("imager-pixel-row-spacing", ps[0], h5); + WriteSingleScalarH5("imager-pixel-col-spacing", ps[1], h5); + } + + WriteOptionalScalarH5("grid-focal-dist-mm", dcm_info.grid_focal_dist_mm, h5); + + WriteOptionalScalarH5("window-center", dcm_info.window_center, h5); + WriteOptionalScalarH5("window-width", dcm_info.window_width, h5); +} + +xreg::DICOMFIleBasicFields xreg::ReadDICOMFieldsH5(const H5::Group& h5) +{ + DICOMFIleBasicFields dcm_info; + + dcm_info.file_path = ReadStringH5("file-path", h5); + + dcm_info.patient_id = ReadStringH5("patient-id", h5); + dcm_info.series_uid = ReadStringH5("series-uid", h5); + dcm_info.study_uid = ReadStringH5("study-uid", h5); + + dcm_info.patient_name = ReadStringH5("patient-name", h5); + + dcm_info.study_time = ReadSingleScalarH5Double("study-time", h5); + dcm_info.series_time = ReadOptionalScalarH5("series-time", h5); + dcm_info.acquisition_time = ReadOptionalScalarH5("acquisition-time", h5); + dcm_info.content_time = ReadOptionalScalarH5("content-time", h5); + + dcm_info.modality = ReadStringH5("modality", h5); + + dcm_info.img_pos_wrt_pat = detail::ReadMatrixH5Helper("img-pos-wrt-pat", h5); + + dcm_info.row_dir = detail::ReadMatrixH5Helper("row-dir", h5); + dcm_info.col_dir = detail::ReadMatrixH5Helper("col-dir", h5); + + dcm_info.row_spacing = ReadSingleScalarH5CoordScalar("row-spacing", h5); + dcm_info.col_spacing = ReadSingleScalarH5CoordScalar("col-spacing", h5); + + dcm_info.num_rows = ReadSingleScalarH5ULong("num-rows", h5); + dcm_info.num_cols = ReadSingleScalarH5ULong("num-cols", h5); + + dcm_info.pat_pos = ReadOptionalStringH5("pat-pos", h5); + + if (ObjectInGroupH5("pat-orient-x", h5) && ObjectInGroupH5("pat-orient-y", h5)) + { + dcm_info.pat_orient = std::array{ ReadStringH5("pat-orient-x", h5), + ReadStringH5("pat-orient-y", h5) }; + } + + dcm_info.study_desc = ReadOptionalStringH5("study-desc", h5); + + dcm_info.series_desc = ReadOptionalStringH5("series-desc", h5); + + dcm_info.image_type = ReadOptionalListOfStrings("image-type", h5); + + dcm_info.manufacturer = ReadStringH5("manufacturer", h5); + + dcm_info.institution_name = ReadOptionalStringH5("institution-name", h5); + + dcm_info.department_name = ReadOptionalStringH5("department-name", h5); + + dcm_info.manufacturers_model_name = ReadOptionalStringH5("manufacturers-model-name", h5); + + dcm_info.sec_cap_dev_manufacturer = ReadOptionalStringH5("sec-cap-dev-manufacturer", h5); + + dcm_info.sec_cap_dev_software_versions = ReadOptionalStringH5("sec-cap-dev-software-versions", h5); + + dcm_info.software_versions = ReadOptionalListOfStrings("software-versions", h5); + + dcm_info.vol_props = ReadOptionalStringH5("vol-props", h5); + + dcm_info.num_frames = ReadOptionalScalarH5("num-frames", h5); + + dcm_info.proto_name = ReadOptionalStringH5("proto-name", h5); + + dcm_info.conv_kernel = ReadOptionalStringH5("conv-kernel", h5); + + dcm_info.body_part_examined = ReadOptionalStringH5("body-part-examined", h5); + + dcm_info.view_position = ReadOptionalStringH5("view-position", h5); + + dcm_info.dist_src_to_det_mm = ReadOptionalScalarH5("dist-src-to-det-mm", h5); + + dcm_info.dist_src_to_pat_mm = ReadOptionalScalarH5("dist-src-to-pat-mm", h5); + + dcm_info.kvp = ReadOptionalScalarH5("kvp", h5); + + dcm_info.tube_current_mA = ReadOptionalScalarH5("tube-current-mA", h5); + + dcm_info.exposure_mAs = ReadOptionalScalarH5("exposure-mAs", h5); + + dcm_info.exposure_muAs = ReadOptionalScalarH5("exposure-muAs", h5); + + dcm_info.exposure_time_ms = ReadOptionalScalarH5("exposure-time-ms", h5); + + dcm_info.dose_area_product_dGy_cm_sq = ReadOptionalScalarH5("dose-area-product-dGy-cm-sq", h5); + + dcm_info.fov_shape = ReadOptionalStringH5("fov-shape", h5); + + if (ObjectInGroupH5("fov-dims", h5)) + { + dcm_info.fov_dims = ReadVectorH5ULong("fov-dims", h5); + } + + if (ObjectInGroupH5("fov-origin-off-rows", h5) && ObjectInGroupH5("fov-origin-off-cols", h5)) + { + dcm_info.fov_origin_off = std::array{ + ReadSingleScalarH5ULong("fov-origin-off-rows", h5), + ReadSingleScalarH5ULong("fov-origin-off-cols", h5) }; + } + + if (ObjectInGroupH5("fov-rot", h5)) + { + dcm_info.fov_rot = static_cast( + ReadSingleScalarH5Int("fov-rot", h5)); + } + + dcm_info.fov_horizontal_flip = ReadOptionalScalarH5("fov-horizontal-flip", h5); + + dcm_info.intensifier_diameter_mm = ReadOptionalScalarH5("intensifier-diameter-mm", h5); + + if (ObjectInGroupH5("imager-pixel-row-spacing", h5) && + ObjectInGroupH5("imager-pixel-col-spacing", h5)) + { + dcm_info.imager_pixel_spacing = std::array{ + ReadSingleScalarH5CoordScalar("imager-pixel-row-spacing", h5), + ReadSingleScalarH5CoordScalar("imager-pixel-col-spacing", h5) }; + } + + dcm_info.grid_focal_dist_mm = ReadOptionalScalarH5("grid-focal-dist-mm", h5); + + dcm_info.window_center = ReadOptionalScalarH5("window-center", h5); + dcm_info.window_width = ReadOptionalScalarH5("window-width", h5); + + return dcm_info; +} + diff --git a/lib/file_formats/xregDICOMUtils.h b/lib/file_formats/xregDICOMUtils.h index 6960b14..a617f1a 100644 --- a/lib/file_formats/xregDICOMUtils.h +++ b/lib/file_formats/xregDICOMUtils.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -31,9 +31,20 @@ #include #include +#include + #include "xregCommon.h" #include "xregFilesystemUtils.h" #include "xregObjWithOStream.h" +#include "xregPerspectiveXform.h" +#include "xregProjData.h" + +// Forward declaration +namespace H5 +{ +class Group; +} +// End Forward declaration namespace xreg { @@ -46,6 +57,13 @@ struct DICOMFIleBasicFields std::string patient_id; std::string series_uid; std::string study_uid; + + std::string patient_name; + + double study_time; + boost::optional series_time; + boost::optional acquisition_time; + boost::optional content_time; std::string modality; @@ -60,47 +78,97 @@ struct DICOMFIleBasicFields unsigned long num_rows; unsigned long num_cols; - bool pat_pos_valid; - std::string pat_pos; + boost::optional pat_pos; + + boost::optional> pat_orient; + + boost::optional study_desc; + + boost::optional series_desc; + + boost::optional> image_type; + + std::string manufacturer; + + boost::optional institution_name; + + boost::optional department_name; + + boost::optional manufacturers_model_name; + + boost::optional sec_cap_dev_manufacturer; + + boost::optional sec_cap_dev_software_versions; - bool pat_orient_valid; - std::array pat_orient; + boost::optional> software_versions; - bool study_desc_valid; - std::string study_desc; + boost::optional vol_props; - bool series_desc_valid; - std::string series_desc; + boost::optional num_frames; - bool image_type_valid; - std::vector image_type; + boost::optional proto_name; - bool sec_cap_dev_manufacturer_valid; - std::string sec_cap_dev_manufacturer; + boost::optional conv_kernel; - bool sec_cap_dev_software_versions_valid; - std::string sec_cap_dev_software_versions; + // Fields that we would like to use from 2D radiographs/fluoro: - bool software_versions_valid; - std::vector software_versions; + boost::optional body_part_examined; - bool vol_props_valid; - std::string vol_props; + boost::optional view_position; - bool num_frames_valid; - unsigned long num_frames; + boost::optional dist_src_to_det_mm; + + boost::optional dist_src_to_pat_mm; + + boost::optional kvp; - bool proto_name_valid; - std::string proto_name; + boost::optional tube_current_mA; - bool conv_kernel_valid; - std::string conv_kernel; + boost::optional exposure_mAs; + + boost::optional exposure_muAs; + + boost::optional exposure_time_ms; + + // units are dGy * cm * cm + boost::optional dose_area_product_dGy_cm_sq; + + boost::optional fov_shape; + + boost::optional> fov_dims; + + boost::optional> fov_origin_off; + + enum FOVRot + { + kZERO = 0, + kNINETY = 90, + kONE_EIGHTY = 180, + kTWO_SEVENTY = 270 + }; + + boost::optional fov_rot; + + boost::optional fov_horizontal_flip; + + boost::optional intensifier_diameter_mm; + + // This is usally populated for 2D X-ray images, e.g. when the standard + // pixel spacing fields are not appropriate as they are required to be + // in "patient space." + // row spacing , col spacing + boost::optional> imager_pixel_spacing; + + boost::optional grid_focal_dist_mm; + + boost::optional window_center; + boost::optional window_width; }; using DICOMFIleBasicFieldsList = std::vector; /// \brief Populates a set of basic DICOM fields from a DICOM file. -void ReadDICOMFileBasicFields(const std::string& dcm_path, DICOMFIleBasicFields* dcm_info); +DICOMFIleBasicFields ReadDICOMFileBasicFields(const std::string& dcm_path); /// \brief Prints a set of basic DICOM fields void PrintDICOMFileBasicFields(const DICOMFIleBasicFields& dcm_info, std::ostream& out, @@ -114,6 +182,10 @@ bool IsVolDICOMFile(const DICOMFIleBasicFields& dcm_info); bool IsMultiFrameDICOMFile(const DICOMFIleBasicFields& dcm_info); +bool IsSecondaryDICOMFile(const DICOMFIleBasicFields& dcm_info); + +bool IsDerivedDICOMFile(const DICOMFIleBasicFields& dcm_info); + /// \brief Stores paths to DICOM files organized by patient ID, study UID, and /// series UID. struct OrganizedDICOMFiles @@ -146,10 +218,15 @@ void GetDICOMFilePathObjsInDir(const std::string& dir, PathList* dcm_paths); /// hierarchy and store in a hierarchy of data structures. /// /// Organized as Patient ID -> Studies for each Patient ID -> Series for each study +/// When the modalities argument is non-empty, then only the specified modalities will +/// be included in the output. void GetOrgainizedDICOMInfos(const std::string& root_dir_path, OrganizedDICOMFiles* org_dcm, const bool inc_localizer = false, - const bool inc_multi_frame_files = false); + const bool inc_multi_frame_files = false, + const bool inc_secondary = true, + const bool inc_derived = true, + const std::vector& modalities = std::vector()); /// \brief Get basic information structs for every DICOM file in a single directory /// @@ -174,10 +251,19 @@ struct ReorderAndCheckDICOMInfos : public ObjWithOStream // if the slice spacing is constant - default 1.0e-6 CoordScalar out_of_plane_spacing_tol = 1.0e-6; + // Perform the reordering and verification. + // The input list of DICOM fields is preserved and the sorted output is populated + // in a separate list via a pointer supplied by the caller. + // This routine returns true when the data is determined to be valid and false otherwise. + // When false is returned, the output sorted list should not be considered to be valid in any way. bool operator()(const DICOMFIleBasicFieldsList& src_infos, DICOMFIleBasicFieldsList* dst_infos); }; +void WriteDICOMFieldsH5(const DICOMFIleBasicFields& dcm_info, H5::Group* h5); + +DICOMFIleBasicFields ReadDICOMFieldsH5(const H5::Group& h5); + } // xreg #endif diff --git a/lib/file_formats/xregH5ProjDataIO.cpp b/lib/file_formats/xregH5ProjDataIO.cpp index b2ff38a..2eddf25 100644 --- a/lib/file_formats/xregH5ProjDataIO.cpp +++ b/lib/file_formats/xregH5ProjDataIO.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -26,10 +26,10 @@ #include +#include "xregDICOMUtils.h" #include "xregHDF5.h" #include "xregHDF5Internal.h" #include "xregH5CamModelIO.h" -#include "xregCIOSFusionDICOM.h" namespace { @@ -90,14 +90,19 @@ void WriteProjDataH5Helper(const std::vector>& proj_data, static_cast(*proj_data[i].rot_to_pat_up), &proj_g); } - if (proj_data[i].orig_meta) + if (proj_data[i].det_spacings_from_orig_meta) { - H5::Group orig_meta_g = proj_g.createGroup("orig-meta"); + WriteSingleScalarH5("det-spacings-from-orig-meta", + *proj_data[i].det_spacings_from_orig_meta, &proj_g); + } + + if (proj_data[i].orig_dcm_meta) + { + H5::Group orig_meta_g = proj_g.createGroup("orig-dcm-meta"); - // For now we are only storing metadata from the CIOS fusion - SetStringAttr("meta-type", "cios-fusion", &orig_meta_g); + SetStringAttr("meta-type", "dicom", &orig_meta_g); - WriteCIOSMetaH5(*proj_data[i].orig_meta, &orig_meta_g); + WriteDICOMFieldsH5(*proj_data[i].orig_dcm_meta, &orig_meta_g); } } } @@ -390,14 +395,23 @@ ReadProjDataHelper(const H5::Group& h5, const bool read_pixels) projs[i].rot_to_pat_up = static_cast( ReadSingleScalarH5Int("rot-to-pat-up", proj_g)); } + + if (ObjectInGroupH5("det-spacings-from-orig-meta", proj_g)) + { + projs[i].det_spacings_from_orig_meta = ReadSingleScalarH5Bool( + "det-spacings-from-orig-meta", proj_g); + } - if (ObjectInGroupH5("orig-meta", proj_g)) + if (ObjectInGroupH5("orig-dcm-meta", proj_g)) { - // TODO: check the "meta-type" attribute after more sensors are added + H5::Group orig_meta_g = proj_g.openGroup("orig-dcm-meta"); + + // double-check that DICOM metadata was written here + xregASSERT(GetStringAttr("meta-type", orig_meta_g) == "dicom"); - projs[i].orig_meta = std::make_shared(); + projs[i].orig_dcm_meta = std::make_shared(); - *projs[i].orig_meta = ReadCIOSMetaH5(proj_g.openGroup("orig-meta")); + *projs[i].orig_dcm_meta = ReadDICOMFieldsH5(orig_meta_g); } } diff --git a/lib/file_formats/xregRadRawProj.cpp b/lib/file_formats/xregRadRawProj.cpp new file mode 100644 index 0000000..59dd3c3 --- /dev/null +++ b/lib/file_formats/xregRadRawProj.cpp @@ -0,0 +1,242 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregRadRawProj.h" + +#include + +#include + +#include "xregAssert.h" +#include "xregFilesystemUtils.h" +#include "xregITKBasicImageUtils.h" +#include "xregRigidUtils.h" +#include "xregStringUtils.h" + +xreg::RadRawProjInfo xreg::ReadRadRawProjInfo(const std::string& file_path) +{ + std::ifstream in(file_path); + + const auto lines = GetNonEmptyLinesFromStream(in); + + const size_type num_lines = lines.size(); + + RadRawProjInfo info; + + bool found_size = false; + bool found_step = false; + bool found_data_type = false; + bool found_tpos = false; + bool found_src_pos = false; + + const char* size_param_str = "Size[pixels]"; + const char* step_param_str = "Step[mm]"; + const char* data_type_param_str = "DataType"; + const char* src_pos_param_str = "SourcePosition[mm]"; + const char* pos_param_str = "TPosition"; + + for (size_type line_idx = 0; line_idx < num_lines; ++line_idx) + { + const auto colon_sep_toks = StringSplit(lines[line_idx], ":"); + + const size_type num_colon_sep_toks = colon_sep_toks.size(); + + if (num_colon_sep_toks > 0) + { + const auto& param_name = colon_sep_toks[0]; + + if (num_colon_sep_toks > 1) + { + const auto& param_val = colon_sep_toks[1]; + + if (param_name == size_param_str) + { + const auto size_toks = StringCast(StringSplit(param_val)); + + xregASSERT(size_toks.size() == 2); + + info.num_cols = size_toks[0]; + info.num_rows = size_toks[1]; + + found_size = true; + } + else if (param_name == step_param_str) + { + const auto step_toks = StringCast(StringSplit(param_val)); + + xregASSERT(step_toks.size() == 2); + + info.col_spacing_mm_per_pixel = step_toks[0]; + info.row_spacing_mm_per_pixel = step_toks[1]; + + found_step = true; + } + else if (param_name == data_type_param_str) + { + info.data_type = StringStrip(param_val); + + found_data_type = true; + } + else if (param_name == src_pos_param_str) + { + const auto pos_toks = StringCast(StringSplit(param_val)); + + xregASSERT(pos_toks.size() == 3); + + info.src_wrt_world(0) = pos_toks[0]; + info.src_wrt_world(1) = pos_toks[1]; + info.src_wrt_world(2) = pos_toks[2]; + + found_src_pos = true; + } + } + else if (param_name == pos_param_str) + { + if ((line_idx + 4) < num_lines) + { + for (size_type r = 0; r < 4; ++r) + { + const auto mat_row_toks = StringCast(StringSplit(lines[line_idx + r + 1])); + xregASSERT(mat_row_toks.size() == 4); + + for (size_type c = 0; c < 4; ++c) + { + info.proj_to_world.matrix()(r,c) = mat_row_toks[c]; + } + } + + line_idx += 4; + + found_tpos = true; + } + } + } + } + + if (!(found_size && found_step && found_data_type && found_tpos && found_src_pos)) + { + std::vector params_not_found; + + if (!found_size) + { + params_not_found.push_back(size_param_str); + } + + if (!found_step) + { + params_not_found.push_back(step_param_str); + } + + if (!found_data_type) + { + params_not_found.push_back(data_type_param_str); + } + + if (!found_tpos) + { + params_not_found.push_back(pos_param_str); + } + + if (!found_src_pos) + { + params_not_found.push_back(src_pos_param_str); + } + + xregThrow("Failed to find required .rad metadata fields: %s", + JoinTokens(params_not_found, " , ").c_str()); + } + + return info; +} + +std::tuple::Pointer> +xreg::ReadRadRawProj(const std::string& rad_file_path) +{ + const auto info = ReadRadRawProjInfo(rad_file_path); + + xregASSERT(info.data_type == "float"); + + auto img = MakeITK2DVol(info.num_cols, info.num_rows); + + { + const std::array tmp_spacing = { info.col_spacing_mm_per_pixel, + info.row_spacing_mm_per_pixel }; + + img->SetSpacing(tmp_spacing.data()); + } + + FileInputStream fin(fmt::format("{}.raw", std::get<0>(Path(rad_file_path).split_ext()))); + + xregASSERT(fin.num_bytes_left() == (sizeof(float) * info.num_cols * info.num_rows)); + + fin.read_remaining_bytes(img->GetBufferPointer()); + + return std::make_tuple(info, img); +} + +xreg::ProjDataF32 xreg::ReadRawProjAsProjData(const std::string& rad_file_path) +{ + ProjDataF32 pd; + + auto rad_data = ReadRadRawProj(rad_file_path); + + const auto& info = std::get<0>(rad_data); + + pd.img = std::get<1>(rad_data); + + // The .rad projective frame has an origin at the (0,0) pixel in the 2D image + // with x axis aligned with increasing columns, y axis with increasing rows, + // and z axis pointing towards the source. + // + // This code moves the origin to be at the source, keeps the axes orientation unchanged, + // and creates an appropriate intrinsic matrix assuming zero shear. + // The extrinsic matrix will be updated to account for translating the origin. + + const FrameTransform rad_world_to_proj = info.proj_to_world.inverse(); + + const Pt3 src_wrt_proj_rad = rad_world_to_proj * info.src_wrt_world; + + const CoordScalar src_to_det_dist = src_wrt_proj_rad(2); + + xregASSERT(src_to_det_dist > 1.0e-6); + + Mat3x3 K = Mat3x3::Identity(); + + K(0,0) = -src_to_det_dist / info.col_spacing_mm_per_pixel; + K(1,1) = -src_to_det_dist / info.row_spacing_mm_per_pixel; + + K(0,2) = src_wrt_proj_rad(0) / info.col_spacing_mm_per_pixel; + K(1,2) = src_wrt_proj_rad(1) / info.row_spacing_mm_per_pixel; + + pd.cam.coord_frame_type = CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z; + + pd.cam.setup(K, TransXYZ4x4(-src_wrt_proj_rad) * rad_world_to_proj.matrix(), + info.num_rows, info.num_cols, + info.col_spacing_mm_per_pixel, info.row_spacing_mm_per_pixel); + + pd.det_spacings_from_orig_meta = true; + + return pd; +} + diff --git a/lib/file_formats/xregRadRawProj.h b/lib/file_formats/xregRadRawProj.h new file mode 100644 index 0000000..9c639e3 --- /dev/null +++ b/lib/file_formats/xregRadRawProj.h @@ -0,0 +1,59 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#ifndef XREGRADRAWPROJ_H_ +#define XREGRADRAWPROJ_H_ + +#include "xregProjData.h" + +namespace xreg +{ + +// Store metadata found in .rad/.raw projection files from the Ljubljana 2D/3D datasets. +// Data is available here: http://lit.fe.uni-lj.si/tools.php?lang=eng +struct RadRawProjInfo +{ + size_type num_cols; + size_type num_rows; + + float col_spacing_mm_per_pixel; + float row_spacing_mm_per_pixel; + + std::string data_type; + + FrameTransform proj_to_world; + + Pt3 src_wrt_world; +}; + +RadRawProjInfo ReadRadRawProjInfo(const std::string& file_path); + +std::tuple::Pointer> ReadRadRawProj(const std::string& rad_file_path); + +ProjDataF32 ReadRawProjAsProjData(const std::string& rad_file_path); + +} // xreg + +#endif + diff --git a/lib/file_formats/xregReadProjDataFromDICOM.cpp b/lib/file_formats/xregReadProjDataFromDICOM.cpp new file mode 100644 index 0000000..78591cb --- /dev/null +++ b/lib/file_formats/xregReadProjDataFromDICOM.cpp @@ -0,0 +1,55 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregReadProjDataFromDICOM.h" +#include "xregReadProjDataFromDICOMDetail.h" + +xreg::ProjDataF32List xreg::ReadProjDataFromDICOMF32(const std::string& dcm_path, + const ReadProjDataFromDICOMParams& params) +{ + return detail::ReadProjDataFromDICOMHelper(dcm_path, params); +} + +xreg::ProjDataF32List xreg::ReadProjDataFromDICOMF32(const std::string& dcm_path, + const std::string& fcsv_path, + const ReadProjDataFromDICOMParams& params) +{ + return detail::ReadProjDataFromDICOMHelper(dcm_path, fcsv_path, params); +} + + +xreg::ProjDataU16List xreg::ReadProjDataFromDICOMU16(const std::string& dcm_path, + const ReadProjDataFromDICOMParams& params) +{ + return detail::ReadProjDataFromDICOMHelper(dcm_path, params); +} + +xreg::ProjDataU16List xreg::ReadProjDataFromDICOMU16(const std::string& dcm_path, + const std::string& fcsv_path, + const ReadProjDataFromDICOMParams& params) +{ + return detail::ReadProjDataFromDICOMHelper(dcm_path, fcsv_path, params); +} + + diff --git a/lib/file_formats/xregReadProjDataFromDICOM.h b/lib/file_formats/xregReadProjDataFromDICOM.h new file mode 100644 index 0000000..b8d8e05 --- /dev/null +++ b/lib/file_formats/xregReadProjDataFromDICOM.h @@ -0,0 +1,89 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#ifndef XREGREADPROJDATAFROMDICOM_H_ +#define XREGREADPROJDATAFROMDICOM_H_ + +#include + +#include "xregProjData.h" + +namespace xreg +{ + +struct ReadProjDataFromDICOMParams +{ + double src_to_det_default = 1000.0; + + double spacing_default = 1.0; + + // attempt to guess the image spacing when the DICOM metadata does not provide an explicit value. + bool guess_spacing = true; + + // will auto-set the projective frame based on modality when no value is provided + boost::optional proj_frame; + + bool no_bayview_check = false; + + // Do not perform any pre-processing to the image pixels - e.g. do NOT flip or rotate the image + // using the DICOM FOV Rotation or FOV Horizontal Flip fields. + bool no_proc = false; + + // This is used for converting landmarks in physical FCSV coordinates to pixel locations. + // When no metadata for row/column spacing is explicitly specified in the DICOM, 3D Slicer + // will typically use a default spacing of 1.0. This field exist in the event that another + // default FCSV spacing needs to be specified. + double fcsv_spacing_default = 1.0; + + // Output stream to print verbose information helpful in debugging, etc. + // A null (e.g. like /dev/null) output stream will be used when nullptr is provided. + std::ostream* vout = nullptr; + + // Output stream to print warnings and error messages. + // std::cerr will be used when given a nullptr. + std::ostream* err_out = nullptr; +}; + +ProjDataF32List ReadProjDataFromDICOMF32(const std::string& dcm_path, + const ReadProjDataFromDICOMParams& params = + ReadProjDataFromDICOMParams()); + +ProjDataF32List ReadProjDataFromDICOMF32(const std::string& dcm_path, + const std::string& fcsv_path, + const ReadProjDataFromDICOMParams& params = + ReadProjDataFromDICOMParams()); + +ProjDataU16List ReadProjDataFromDICOMU16(const std::string& dcm_path, + const ReadProjDataFromDICOMParams& params = + ReadProjDataFromDICOMParams()); + +ProjDataU16List ReadProjDataFromDICOMU16(const std::string& dcm_path, + const std::string& fcsv_path, + const ReadProjDataFromDICOMParams& params = + ReadProjDataFromDICOMParams()); + +} // xreg + +#endif + diff --git a/lib/file_formats/xregReadProjDataFromDICOMDetail.h b/lib/file_formats/xregReadProjDataFromDICOMDetail.h new file mode 100644 index 0000000..90140a4 --- /dev/null +++ b/lib/file_formats/xregReadProjDataFromDICOMDetail.h @@ -0,0 +1,519 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#ifndef XREGREADPROJDATAFROMDICOMDETAIL_H_ +#define XREGREADPROJDATAFROMDICOMDETAIL_H_ + +#include +#include + +#include + +#include "xregAnatCoordFrames.h" +#include "xregCIOSFusionDICOM.h" +#include "xregDICOMUtils.h" +#include "xregFCSVUtils.h" +#include "xregITKIOUtils.h" +#include "xregITKOpenCVUtils.h" +#include "xregLandmarkMapUtils.h" +#include "xregOpenCVUtils.h" +#include "xregReadProjDataFromDICOM.h" + +namespace xreg +{ +namespace detail +{ + +template +std::vector> +ReadProjDataFromDICOMHelper(const std::string& dcm_path, const ReadProjDataFromDICOMParams& params) +{ + boost::iostreams::stream null_ostream((boost::iostreams::null_sink())); + + std::ostream& vout = params.vout ? *params.vout : null_ostream; + std::ostream& err_out = params.err_out ? *params.err_out : std::cerr; + + vout << "reading DICOM metadata..." << std::endl; + const auto dcm_info = ReadDICOMFileBasicFields(dcm_path); + + vout << " input modality: " << dcm_info.modality << std::endl; + + const bool modality_is_xa = dcm_info.modality == "XA"; + const bool modality_is_dx = dcm_info.modality == "DX"; + const bool modality_is_cr = dcm_info.modality == "CR"; + const bool modality_is_rf = dcm_info.modality == "RF"; + + if (!(modality_is_xa || modality_is_dx || modality_is_cr || modality_is_rf)) + { + err_out << "WARNING: unexpected modality: " << dcm_info.modality << std::endl; + } + + CameraModel::CameraCoordFrame proj_frame; + + if (!params.proj_frame) + { + vout << "automatically selecting proj. frame Z direction using modality..." << std::endl; + + proj_frame = (modality_is_xa || modality_is_rf) ? + CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z : + CameraModel::kORIGIN_AT_FOCAL_PT_DET_POS_Z; + } + else + { + vout << "using specified value for camera coord frame: " + << static_cast(*params.proj_frame) << std::endl; + + proj_frame = *params.proj_frame; + } + + vout << "setting up camera model..." << std::endl; + + float src_to_det_to_use = static_cast(params.src_to_det_default); + + float row_spacing_to_use = static_cast(params.spacing_default); + float col_spacing_to_use = row_spacing_to_use; + + if (dcm_info.dist_src_to_det_mm) + { + src_to_det_to_use = *dcm_info.dist_src_to_det_mm; + } + else + { + err_out << "WARNING: source to detector field not present in DICOM, will use default value of " + << params.src_to_det_default << std::endl; + } + + bool spacing_in_meta = true; + + // prefer to use the imager spacing field when available + if (dcm_info.imager_pixel_spacing) + { + vout << "using imager pixel spacing field" << std::endl; + + const auto& s = *dcm_info.imager_pixel_spacing; + + row_spacing_to_use = s[0]; + col_spacing_to_use = s[1]; + } + else if ((dcm_info.row_spacing > 1.0e-6) && (dcm_info.col_spacing > 1.0e-6)) + { + // next, use the image pixel spacing field - this is less preferred than the + // imager spacing as this field is supposed to be defined with respect to a + // patient coordinate frame, which does not make sense for a 2D radiograph + + vout << "using image pixel spacing..." << std::endl; + + row_spacing_to_use = dcm_info.row_spacing; + col_spacing_to_use = dcm_info.col_spacing; + } + else + { + // No other tags explicitly specify the spacing, try and guess using some other metadata fields. + + spacing_in_meta = false; + + bool guess_made = false; + + if (params.guess_spacing) + { + vout << "pixel spacing metadata not available - attempting to guess..." << std::endl; + + if (dcm_info.fov_shape) + { + const auto& fov_shape_str = *dcm_info.fov_shape; + + vout << " FOV shape available: " << fov_shape_str << std::endl; + + if (dcm_info.fov_dims) + { + const auto& fov_dims = *dcm_info.fov_dims; + + const size_type num_fov_dims = fov_dims.size(); + + unsigned long num_rows_for_guess = dcm_info.num_rows; + unsigned long num_cols_for_guess = dcm_info.num_cols; + + if (dcm_info.fov_origin_off) + { + const auto& fov_origin_off = *dcm_info.fov_origin_off; + + vout << "FOV origin offset available: [ " << fov_origin_off[0] << " , " + << fov_origin_off[1] << " ]" << std::endl; + + num_rows_for_guess -= 2 * fov_origin_off[0]; + num_cols_for_guess -= 2 * fov_origin_off[1]; + + vout << " number of rows used for spacing guess: " << num_rows_for_guess << std::endl; + vout << " number of cols used for spacing guess: " << num_cols_for_guess << std::endl; + } + + if (fov_shape_str == "ROUND") + { + if (num_fov_dims == 1) + { + if (dcm_info.num_cols != dcm_info.num_rows) + { + err_out << "WARNING: number of image rows and columns are not equal!" + "Guessed pixels spacings may have substantial errors!" << std::endl; + } + + row_spacing_to_use = static_cast(fov_dims[0]) / + std::max(num_cols_for_guess, num_rows_for_guess); + col_spacing_to_use = row_spacing_to_use; + + vout << " using round FOV diameter of " << fov_dims[0] + << " mm to guess isotropic spacing of " + << fmt::format("{:.3f}", row_spacing_to_use) << " mm/pixel" << std::endl; + + guess_made = true; + } + else + { + err_out << "expected ROUND FOV dims to have length 1, got: " << num_fov_dims << std::endl; + } + } + else if (fov_shape_str == "RECTANGLE") + { + if (num_fov_dims == 2) + { + col_spacing_to_use = static_cast(fov_dims[0]) / num_cols_for_guess; + row_spacing_to_use = static_cast(fov_dims[1]) / num_rows_for_guess; + + vout << " using rect FOV row length of " << fov_dims[0] + << " mm to guess column spacing of " + << fmt::format("{:.3f}", col_spacing_to_use) << " mm/pixel" << std::endl; + + vout << " using rect FOV column length of " << fov_dims[0] + << " mm to guess row spacing of " + << fmt::format("{:.3f}", row_spacing_to_use) << " mm/pixel" << std::endl; + + guess_made = true; + } + else + { + err_out << "expected RECTANGLE FOV dims to have length 2, got: " + << num_fov_dims << std::endl; + } + } + else + { + err_out << "unsupported/unknown FOV shape string: " << fov_shape_str << std::endl; + } + } + else + { + vout << " FOV dims not available" << std::endl; + } + } + + if (!guess_made) + { + if (dcm_info.intensifier_diameter_mm) + { + if (dcm_info.num_cols != dcm_info.num_rows) + { + err_out << "WARNING: number of image rows and columns are not equal!" + "Guessed pixels spacings may have substantial errors!" << std::endl; + } + + const float d = static_cast(*dcm_info.intensifier_diameter_mm); + + row_spacing_to_use = d / std::max(dcm_info.num_cols, dcm_info.num_rows); + col_spacing_to_use = row_spacing_to_use; + + vout << "using intensifier diameter of " << d + << " mm to guess isotropic pixel spacing of " << row_spacing_to_use + << " mm / pixel" << std::endl; + + guess_made = true; + } + } + } + + if (!guess_made) + { + vout << "spacing not found in metadata, using default spacing: " + << params.spacing_default << std::endl; + } + } + + CameraModel cam; + + cam.coord_frame_type = proj_frame; + + const bool is_bayview_cios_dcm = (dcm_info.manufacturer == "SIEMENS") && + (dcm_info.institution_name && + (*dcm_info.institution_name == "Johns Hopkins Univ")) && + (dcm_info.department_name && + (*dcm_info.department_name == "Applied Physics Lab")) && + (dcm_info.manufacturers_model_name && + (*dcm_info.manufacturers_model_name == "Fluorospot Compact S1")) && + dcm_info.dist_src_to_det_mm; + + if (params.no_bayview_check || !is_bayview_cios_dcm) + { + vout << "setting camera model with naive intrinsics and identity extrinsics..." << std::endl; + + cam.setup(src_to_det_to_use, + dcm_info.num_rows, dcm_info.num_cols, + row_spacing_to_use, col_spacing_to_use); + } + else + { + vout << "bayview file detected, setting up camera model with calculated extrinsics..." << std::endl; + + if (cam.coord_frame_type != CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z) + { + err_out << "WARNING: C-arm projective frame type does not match the expected value!\n " + " Expected: " << static_cast(CameraModel::kORIGIN_AT_FOCAL_PT_DET_NEG_Z) + << ", have: " << static_cast(cam.coord_frame_type) + << std::endl; + } + + const Mat3x3 intrins = MakeNaiveIntrins(*dcm_info.dist_src_to_det_mm, + dcm_info.num_rows, dcm_info.num_cols, + row_spacing_to_use, + col_spacing_to_use, + true); + + cam.setup(intrins, CIOSFusionCBCTExtrins(), + dcm_info.num_rows, dcm_info.num_cols, + row_spacing_to_use, col_spacing_to_use); + } + + const size_type num_frames = dcm_info.num_frames ? * dcm_info.num_frames : 1; + + std::vector> pd(num_frames); + + if (num_frames == 1) + { + vout << "1 frame - reading 2D image pixel data from DICOM..." << std::endl; + pd[0].img = ReadDICOM2DFromDisk(dcm_path); + } + else + { + vout << num_frames << " frames - reading 3D image pixel data from DICOM..." << std::endl; + auto frames = ReadDICOM3DFromDisk(dcm_path); + + const auto spacing_vol = frames->GetSpacing(); + + const std::array spacing_slice = { spacing_vol[0], spacing_vol[1] }; + + const auto origin_vol = frames->GetOrigin(); + + const std::array origin_slice = { origin_vol[0], origin_vol[1] }; + + vout << " converting in-plane slices to individual projection frames..." << std::endl; + + const auto* cur_frame_buf = frames->GetBufferPointer(); + + const size_type num_pix_per_frame = cam.num_det_cols * cam.num_det_rows; + + for (size_type i = 0; i < num_frames; ++i, cur_frame_buf += num_pix_per_frame) + { + auto dst_frame = MakeITK2DVol(cam.num_det_cols, cam.num_det_rows); + + dst_frame->SetSpacing(spacing_slice.data()); + dst_frame->SetOrigin(origin_slice.data()); + + std::copy(cur_frame_buf, cur_frame_buf + num_pix_per_frame, dst_frame->GetBufferPointer()); + + pd[i].img = dst_frame; + } + } + + { + auto img_spacing = pd[0].img->GetSpacing(); + + if (std::abs(img_spacing[0] - cam.det_col_spacing) > 1.0e-3) + { + err_out << "WARNING: Image column spacing (" << img_spacing[0] + <<") differs from camera model column spacings (" + << cam.det_col_spacing + << "). Image values will be updated to match camera model." + << std::endl; + } + + if (std::abs(img_spacing[1] - cam.det_row_spacing) > 1.0e-3) + { + err_out << "WARNING: Image row spacing (" << img_spacing[1] + <<") differs from camera model row spacings (" + << cam.det_row_spacing + << "). Image values will be updated to match camera model." + << std::endl; + } + } + + // Always prefer the spacing obtained by interpreting DICOM fields + const std::array spacing_to_use = { cam.det_col_spacing, cam.det_row_spacing }; + + auto orig_dcm_meta = std::make_shared(); + *orig_dcm_meta = dcm_info; + + const auto dcm_rot = dcm_info.fov_rot ? *dcm_info.fov_rot : DICOMFIleBasicFields::kZERO; + + const bool do_horiz_flip = dcm_info.fov_horizontal_flip && (*dcm_info.fov_horizontal_flip); + + for (size_type i = 0; i < num_frames; ++i) + { + pd[i].cam = cam; + + // Always prefer the spacing obtained by interpreting DICOM fields + pd[i].img->SetSpacing(spacing_to_use.data()); + + pd[i].det_spacings_from_orig_meta = spacing_in_meta; + + pd[i].orig_dcm_meta = orig_dcm_meta; + + if (!params.no_proc) + { + cv::Mat img_ocv = ShallowCopyItkToOpenCV(pd[i].img.GetPointer()); + + if (dcm_rot != DICOMFIleBasicFields::kZERO) + { + if (dcm_rot == DICOMFIleBasicFields::kNINETY) + { + xregASSERT(pd[i].cam.num_det_rows == pd[i].cam.num_det_cols); + + cv::Mat tmp = img_ocv.clone(); + cv::transpose(tmp, img_ocv); + FlipImageColumns(&img_ocv); + } + else if (dcm_rot == DICOMFIleBasicFields::kONE_EIGHTY) + { + FlipImageRows(&img_ocv); + FlipImageColumns(&img_ocv); + } + else if (dcm_rot == DICOMFIleBasicFields::kTWO_SEVENTY) + { + xregASSERT(pd[i].cam.num_det_rows == pd[i].cam.num_det_cols); + + cv::Mat tmp = img_ocv.clone(); + cv::transpose(tmp, img_ocv); + FlipImageRows(&img_ocv); + } + } + + if (do_horiz_flip) + { + FlipImageColumns(&img_ocv); + } + } + } + + return pd; +} + +template +std::vector> +ReadProjDataFromDICOMHelper(const std::string& dcm_path, const std::string& fcsv_path, + const ReadProjDataFromDICOMParams& params) +{ + boost::iostreams::stream null_ostream((boost::iostreams::null_sink())); + + std::ostream& vout = params.vout ? *params.vout : null_ostream; + + auto pd = ReadProjDataFromDICOMHelper(dcm_path, params); + + if (!fcsv_path.empty()) + { + vout << "reading landmarks from FCSV and converting to pixels..." << std::endl; + auto lands_3d = ReadFCSVFileNamePtMap(fcsv_path); + + ConvertRASToLPS(&lands_3d); + + xregASSERT(!pd.empty()); + + // We need to use the pixel spacing that was used by 3D Slicer to create the FCSV file + // when converting landmarks to pixel indices. The FCSV pixel spacing will be equal to + // the spacing used in our projection data camera model when the spacing was explicitly + // available in the DICOM metadata. Otherwise, we need to have the FCSV spacing + // specified. The FCSV spacing is typically 1.0 when no metadata is provided in the + // DICOM. + + const bool det_spacings_from_orig_meta = *pd[0].det_spacings_from_orig_meta; + + const auto lands = PhysPtsToInds(DropPtDim(lands_3d, 2), + det_spacings_from_orig_meta ? + pd[0].cam.det_col_spacing : + static_cast(params.fcsv_spacing_default), + det_spacings_from_orig_meta ? + pd[0].cam.det_row_spacing : + static_cast(params.fcsv_spacing_default)); + + for (auto& p : pd) + { + p.landmarks = lands; + + if (!params.no_proc && p.orig_dcm_meta) + { + const auto& dcm_info = *p.orig_dcm_meta; + + const auto dcm_rot = dcm_info.fov_rot ? *dcm_info.fov_rot : DICOMFIleBasicFields::kZERO; + + const bool do_horiz_flip = dcm_info.fov_horizontal_flip && (*dcm_info.fov_horizontal_flip); + + for (auto& lkv: p.landmarks) + { + auto& l = lkv.second; + + if (dcm_rot != DICOMFIleBasicFields::kZERO) + { + if (dcm_rot == DICOMFIleBasicFields::kNINETY) + { + std::swap(l(0), l(1)); + } + else if (dcm_rot == DICOMFIleBasicFields::kONE_EIGHTY) + { + l(0) = p.cam.num_det_cols - 1 - l(0); + l(1) = p.cam.num_det_rows - 1 - l(1); + } + else if (dcm_rot == DICOMFIleBasicFields::kTWO_SEVENTY) + { + std::swap(l(0), l(1)); + l(1) = p.cam.num_det_rows - 1 - l(1); + } + } + + if (do_horiz_flip) + { + l(0) = p.cam.num_det_cols - 1 - l(0); + } + } + } + } + } + else + { + vout << "empty path provided for FCSV file - no landmarks will be populated" << std::endl; + } + + return pd; +} + +} // detail +} // xreg + +#endif + diff --git a/lib/file_formats/xregStaRawVol.cpp b/lib/file_formats/xregStaRawVol.cpp new file mode 100644 index 0000000..1e0fe83 --- /dev/null +++ b/lib/file_formats/xregStaRawVol.cpp @@ -0,0 +1,164 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include "xregStaRawVol.h" + +#include + +#include + +#include "xregAssert.h" +#include "xregFilesystemUtils.h" +#include "xregITKBasicImageUtils.h" +#include "xregStringUtils.h" + +xreg::StaVolInfo xreg::ReadStaVolInfo(const std::string& sta_file_path) +{ + std::ifstream in(sta_file_path); + + const auto lines = GetNonEmptyLinesFromStream(in); + + const size_type num_lines = lines.size(); + + StaVolInfo info; + + bool found_size = false; + bool found_step = false; + bool found_data_type = false; + bool found_tpos = false; + + for (size_type line_idx = 0; line_idx < num_lines; ++line_idx) + { + const auto colon_sep_toks = StringSplit(lines[line_idx], ":"); + + const size_type num_colon_sep_toks = colon_sep_toks.size(); + + if (num_colon_sep_toks > 0) + { + const auto& param_name = colon_sep_toks[0]; + + if (num_colon_sep_toks > 1) + { + const auto& param_val = colon_sep_toks[1]; + + if (param_name == "Size[pixels]") + { + const auto size_toks = StringCast(StringSplit(param_val)); + + xregASSERT(size_toks.size() == 3); + + info.dims[0] = size_toks[0]; + info.dims[1] = size_toks[1]; + info.dims[2] = size_toks[2]; + + found_size = true; + } + else if (param_name == "Step[mm]") + { + const auto step_toks = StringCast(StringSplit(param_val)); + + xregASSERT(step_toks.size() == 3); + + info.spacing[0] = step_toks[0]; + info.spacing[1] = step_toks[1]; + info.spacing[2] = step_toks[2]; + + found_step = true; + } + else if (param_name == "DataType") + { + info.data_type_str = StringStrip(param_val); + + found_data_type = true; + } + } + else if (param_name == "TPosition") + { + if ((line_idx + 4) < num_lines) + { + for (size_type r = 0; r < 4; ++r) + { + const auto mat_row_toks = StringCast(StringSplit(lines[line_idx + r + 1])); + xregASSERT(mat_row_toks.size() == 4); + + for (size_type c = 0; c < 4; ++c) + { + info.vol_to_world.matrix()(r,c) = mat_row_toks[c]; + } + } + + line_idx += 4; + + found_tpos = true; + } + } + } + } + + if (!(found_size && found_step && found_data_type && found_tpos)) + { + xregThrow("Failed to find required .sta metadata fields!"); + } + + return info; +} + +itk::Image::Pointer xreg::ReadStaRawVol(const std::string& sta_file_path) +{ + // read in the metadata from the STA text file + const auto sta_info = ReadStaVolInfo(sta_file_path); + + // only supporting uint8 for now as all of the data is of this type + xregASSERT(sta_info.data_type_str == "uint8"); + + // allocate the ITK image object that will be returned, this allocates the pixel buffer + auto vol = MakeITK3DVol(sta_info.dims[0], sta_info.dims[1], sta_info.dims[2]); + + // set the important metadata in the ITK image + + vol->SetSpacing(sta_info.spacing.data()); + + SetITKOriginPoint(vol.GetPointer(), Pt3(sta_info.vol_to_world.matrix().block(0,3,3,1))); + + SetITKDirectionMatrix(vol.GetPointer(), Mat3x3(sta_info.vol_to_world.matrix().block(0,0,3,3))); + + // read in the pixels from the .raw file, casting them from uint8 to float + { + FileInputStream fin(fmt::format("{}.raw", std::get<0>(Path(sta_file_path).split_ext()))); + + const size_type tot_num_voxels = sta_info.dims[0] * sta_info.dims[1] * sta_info.dims[2]; + + xregASSERT(fin.num_bytes_left() == tot_num_voxels); + + std::vector vol_buf_as_uint8(tot_num_voxels); + + fin.read_remaining_bytes(vol_buf_as_uint8.data()); + + std::transform(vol_buf_as_uint8.begin(), vol_buf_as_uint8.end(), vol->GetBufferPointer(), + [] (const FileInputStream::uint8& x) { return static_cast(x); }); + } + + return vol; +} + diff --git a/lib/file_formats/xregStaRawVol.h b/lib/file_formats/xregStaRawVol.h new file mode 100644 index 0000000..01ec3dd --- /dev/null +++ b/lib/file_formats/xregStaRawVol.h @@ -0,0 +1,53 @@ +/* + * MIT License + * + * Copyright (c) 2021 Robert Grupp + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#ifndef XREGSTARAWVOL_H_ +#define XREGSTARAWVOL_H_ + +#include "xregCommon.h" + +namespace xreg +{ + +// Store metadata found in .sta/.raw volume files from the Ljubljana 2D/3D datasets. +// Data is available here: http://lit.fe.uni-lj.si/tools.php?lang=eng +struct StaVolInfo +{ + std::array dims; + + std::array spacing; + + std::string data_type_str; + + FrameTransform vol_to_world; +}; + +StaVolInfo ReadStaVolInfo(const std::string& sta_file_path); + +itk::Image::Pointer ReadStaRawVol(const std::string& sta_file_path); + +} // xreg + +#endif + diff --git a/lib/file_formats/xregWriteVideo.cpp b/lib/file_formats/xregWriteVideo.cpp index 4d8485f..fba8ada 100644 --- a/lib/file_formats/xregWriteVideo.cpp +++ b/lib/file_formats/xregWriteVideo.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020,2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -39,7 +39,11 @@ #include "xregAssert.h" #include "xregFilesystemUtils.h" - + +#ifdef __APPLE__ +#include "xregAppleAVFoundation.h" +#endif + void xreg::WriteImageFramesToVideo::write(const std::vector& frames) { for (const cv::Mat& f : frames) @@ -102,7 +106,7 @@ void xreg::WriteImageFramesToVideoWithFFMPEG::open() dst_vid_path, // output file bp::std_in < p->ffmpeg_p_in, bp::std_out > bp::null, - bp::std_err > bp::null)); + bp::std_err > bp::null)); } void xreg::WriteImageFramesToVideoWithFFMPEG::close() @@ -177,11 +181,19 @@ std::unique_ptr xreg::GetWriteImageFramesToVideo( else #endif { +#ifdef __APPLE__ + std::cerr << "WARNING: could not find FFMPEG executable, falling back to " + "Apple AVFoundation video writer!" << std::endl; + writer.reset(new WriteImageFramesToVideoAppleAVF); +#else + #ifndef _WIN32 - std::cerr << "WARNING: could not find FFMPEG executable, falling back to OpenCV video writer!" << std::endl; + std::cerr << "WARNING: could not find FFMPEG executable, falling back to " + "OpenCV video writer!" << std::endl; #endif writer.reset(new WriteImageFramesToVideoWithOpenCV); +#endif } return writer; @@ -189,17 +201,70 @@ std::unique_ptr xreg::GetWriteImageFramesToVideo( void xreg::WriteAllImageFramesToVideo(const std::string& vid_path, const std::vector& frames, - const double fps) + const double fps_or_len, + const bool is_fps) +{ + if (!frames.empty()) + { + auto writer = GetWriteImageFramesToVideo(); + + writer->dst_vid_path = vid_path; + writer->fps = is_fps ? fps_or_len : (frames.size() / fps_or_len); + + writer->open(); + + writer->write(frames); + + writer->close(); + } + else + { + xregThrow("No frames provided to create video!"); + } +} + +void xreg::WriteImageFilesToVideo(const std::string& vid_path, + const std::vector& img_paths, + const double fps_or_len, + const bool is_fps) { - auto writer = GetWriteImageFramesToVideo(); + std::vector frames; - writer->dst_vid_path = vid_path; - writer->fps = fps; + const size_type num_frames = img_paths.size(); + + frames.reserve(num_frames); + + for (size_type i = 0; i < num_frames; ++i) + { + frames.push_back(cv::imread(img_paths[i])); + } + + WriteAllImageFramesToVideo(vid_path, frames, fps_or_len, is_fps); +} + +void xreg::WriteDirOfImagesToVideo(const std::string& vid_path, + const std::string& img_dir, + const bool lex_sort, + const std::vector& img_exts, + const double fps_or_len, + const bool is_fps) +{ + FileExtensions file_exts; - writer->open(); + for (const auto& ext : img_exts) + { + file_exts.add(ext); + } + + std::vector img_paths; - writer->write(frames); - - writer->close(); + GetFilePathsFromDir(img_dir, &img_paths, file_exts); + + if (lex_sort) + { + std::sort(img_paths.begin(), img_paths.end()); + } + + WriteImageFilesToVideo(vid_path, img_paths, fps_or_len, is_fps); } diff --git a/lib/file_formats/xregWriteVideo.h b/lib/file_formats/xregWriteVideo.h index 6631a60..8f4ae9b 100644 --- a/lib/file_formats/xregWriteVideo.h +++ b/lib/file_formats/xregWriteVideo.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020,2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -118,9 +118,37 @@ class WriteImageFramesToVideoWithOpenCV : public WriteImageFramesToVideo std::unique_ptr GetWriteImageFramesToVideo(); +// The final two arguments are used to determine the speed or length of the video. +// When is_fps == true, then fps_or_len represents the desired frames per second +// of the output video. +// When is_fps == false, then fps_or_len represents the desired length of the output +// video in seconds. void WriteAllImageFramesToVideo(const std::string& vid_path, const std::vector& frames, - const double fps = 10.0); + const double fps_or_len = 10.0, + const bool is_fps = true); + +// The final two arguments are used to determine the speed or length of the video. +// When is_fps == true, then fps_or_len represents the desired frames per second +// of the output video. +// When is_fps == false, then fps_or_len represents the desired length of the output +// video in seconds. +void WriteImageFilesToVideo(const std::string& vid_path, + const std::vector& img_paths, + const double fps_or_len = 10.0, + const bool is_fps = true); + +// The final two arguments are used to determine the speed or length of the video. +// When is_fps == true, then fps_or_len represents the desired frames per second +// of the output video. +// When is_fps == false, then fps_or_len represents the desired length of the output +// video in seconds. +void WriteDirOfImagesToVideo(const std::string& vid_path, + const std::string& img_dir, + const bool lex_sort = false, + const std::vector& img_exts = { ".png" }, + const double fps_or_len = 10.0, + const bool is_fps = true); } // xreg diff --git a/lib/hdf5/xregHDF5.cpp b/lib/hdf5/xregHDF5.cpp index 6ca8c8a..9e939ef 100644 --- a/lib/hdf5/xregHDF5.cpp +++ b/lib/hdf5/xregHDF5.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -27,6 +27,25 @@ #include "xregHDF5Internal.h" #include "xregITKBasicImageUtils.h" +void xreg::SetScalarAttr(const std::string& key, const long val, H5::Group* h5) +{ + const auto dt = LookupH5DataType(); + + H5::Attribute attr = h5->createAttribute(key, dt, H5S_SCALAR); + + attr.write(dt, &val); +} + +long xreg::GetScalarLongAttr(const std::string& key, const H5::Group& h5) +{ + const H5::Attribute attr = h5.openAttribute(key); + + long val; + attr.read(attr.getDataType(), &val); + + return val; +} + H5::DataType xreg::GetH5StringDataType() { return LookupH5DataType(); @@ -34,7 +53,7 @@ H5::DataType xreg::GetH5StringDataType() H5::DataType xreg::GetH5StringDataType(const std::string& s) { - return H5::StrType(H5::PredType::C_S1, s.size()); + return H5::StrType(H5::PredType::C_S1, std::max(std::string::size_type(1), s.size())); } bool xreg::SetStringAttr(const std::string& key, const std::string& val, H5::Group* h5) diff --git a/lib/hdf5/xregHDF5.h b/lib/hdf5/xregHDF5.h index 4af20cb..66a91c8 100644 --- a/lib/hdf5/xregHDF5.h +++ b/lib/hdf5/xregHDF5.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -112,7 +112,10 @@ inline H5::DataType LookupH5DataType() template <> inline H5::DataType LookupH5DataType() { - return H5::StrType(H5::PredType::C_S1); + H5::StrType str_type(H5::PredType::C_S1); + str_type.setSize(H5T_VARIABLE); + + return str_type; } #ifdef _WIN32 @@ -125,6 +128,10 @@ inline H5::DataType LookupH5DataType() #endif +void SetScalarAttr(const std::string& key, const long val, H5::Group* h5); + +long GetScalarLongAttr(const std::string& key, const H5::Group& h5); + H5::DataType GetH5StringDataType(); H5::DataType GetH5StringDataType(const std::string& s); diff --git a/lib/image/xregProjData.cpp b/lib/image/xregProjData.cpp index 831bd81..c920bc5 100644 --- a/lib/image/xregProjData.cpp +++ b/lib/image/xregProjData.cpp @@ -37,16 +37,53 @@ using namespace xreg; template ProjData -DownsampleProjDataHelper(const ProjData& src_proj, const CoordScalar ds_factor) +DownsampleProjDataHelper(const ProjData& src_proj, const CoordScalar ds_factor, + const bool force_even_dims) { ProjData dst_proj; - dst_proj.cam = DownsampleCameraModel(src_proj.cam, ds_factor); + dst_proj.cam = DownsampleCameraModel(src_proj.cam, ds_factor, force_even_dims); if (src_proj.img) { dst_proj.img = DownsampleImage(src_proj.img.GetPointer(), ds_factor); + if (force_even_dims) + { + using ProjImg = typename ProjData::Proj; + using ROIFilter = itk::RegionOfInterestImageFilter; + + auto ds_sz = dst_proj.img->GetLargestPossibleRegion().GetSize(); + + const bool cols_are_odd = ds_sz[0] % 2; + const bool rows_are_odd = ds_sz[1] % 2; + + if (cols_are_odd || rows_are_odd) + { + if (cols_are_odd) + { + --ds_sz[0]; + } + + if (rows_are_odd) + { + --ds_sz[1]; + } + + typename ProjImg::RegionType roi; + roi.SetSize(ds_sz); + roi.SetIndex(0,0); + roi.SetIndex(1,0); + + auto roi_filter = ROIFilter::New(); + roi_filter->SetInput(dst_proj.img); + roi_filter->SetRegionOfInterest(roi); + roi_filter->Update(); + + dst_proj.img = roi_filter->GetOutput(); + } + } + const auto ds_sz = dst_proj.img->GetLargestPossibleRegion().GetSize(); xregASSERT(ds_sz[0] == dst_proj.cam.num_det_cols); xregASSERT(ds_sz[1] == dst_proj.cam.num_det_rows); @@ -65,14 +102,15 @@ DownsampleProjDataHelper(const ProjData& src_proj, const CoordScal template std::vector> -DownsampleProjDataHelper(const std::vector>& src_projs, const CoordScalar ds_factor) +DownsampleProjDataHelper(const std::vector>& src_projs, const CoordScalar ds_factor, + const bool force_even_dims) { std::vector> dst_projs; dst_projs.reserve(src_projs.size()); for (const auto& src_proj : src_projs) { - dst_projs.push_back(DownsampleProjDataHelper(src_proj, ds_factor)); + dst_projs.push_back(DownsampleProjDataHelper(src_proj, ds_factor, force_even_dims)); } return dst_projs; @@ -277,37 +315,43 @@ MakeImageFromCam(const CameraModel& cam) } // un-named -xreg::ProjDataF32 xreg::DownsampleProjData(const ProjDataF32& src_proj, const CoordScalar ds_factor) +xreg::ProjDataF32 xreg::DownsampleProjData(const ProjDataF32& src_proj, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_proj, ds_factor); + return DownsampleProjDataHelper(src_proj, ds_factor, force_even_dims); } -xreg::ProjDataU16 xreg::DownsampleProjData(const ProjDataU16& src_proj, const CoordScalar ds_factor) +xreg::ProjDataU16 xreg::DownsampleProjData(const ProjDataU16& src_proj, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_proj, ds_factor); + return DownsampleProjDataHelper(src_proj, ds_factor, force_even_dims); } -xreg::ProjDataU8 xreg::DownsampleProjData(const ProjDataU8& src_proj, const CoordScalar ds_factor) +xreg::ProjDataU8 xreg::DownsampleProjData(const ProjDataU8& src_proj, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_proj, ds_factor); + return DownsampleProjDataHelper(src_proj, ds_factor, force_even_dims); } xreg::ProjDataF32List -xreg::DownsampleProjData(const ProjDataF32List& src_projs, const CoordScalar ds_factor) +xreg::DownsampleProjData(const ProjDataF32List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_projs, ds_factor); + return DownsampleProjDataHelper(src_projs, ds_factor, force_even_dims); } xreg::ProjDataU16List -xreg::DownsampleProjData(const ProjDataU16List& src_projs, const CoordScalar ds_factor) +xreg::DownsampleProjData(const ProjDataU16List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_projs, ds_factor); + return DownsampleProjDataHelper(src_projs, ds_factor, force_even_dims); } xreg::ProjDataU8List -xreg::DownsampleProjData(const ProjDataU8List& src_projs, const CoordScalar ds_factor) +xreg::DownsampleProjData(const ProjDataU8List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims) { - return DownsampleProjDataHelper(src_projs, ds_factor); + return DownsampleProjDataHelper(src_projs, ds_factor, force_even_dims); } std::vector diff --git a/lib/image/xregProjData.h b/lib/image/xregProjData.h index 0bb9b84..f7e0ee3 100644 --- a/lib/image/xregProjData.h +++ b/lib/image/xregProjData.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020-2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -35,8 +35,8 @@ namespace xreg { // Forward declarations: -struct CIOSFusionDICOMInfo; - +struct DICOMFIleBasicFields; + enum class ProjDataRotToPatUp { kZERO = 0, @@ -68,10 +68,17 @@ struct ProjData // of the image and the inferior portion is approximately located in // the bottom of the image. boost::optional rot_to_pat_up; - - // Original metadata from the sensor for this image - does not need + + // Indicates that the detector spacings specified in cam were explicitly + // defined from metadata fields in the original source. Examples where this + // can be false are when converting from DICOM and the user overrides the + // spacing value manually or the spacing values are guessed from other + // metadata values (e.g. the detector diameter). + boost::optional det_spacings_from_orig_meta; + + // Original DICOM metadata this image - does not need // to be set, e.g. for the case of simulated data - std::shared_ptr orig_meta; + std::shared_ptr orig_dcm_meta; }; using ProjDataF32 = ProjData; @@ -82,13 +89,19 @@ using ProjDataF32List = std::vector; using ProjDataU16List = std::vector; using ProjDataU8List = std::vector; -ProjDataF32 DownsampleProjData(const ProjDataF32& src_proj, const CoordScalar ds_factor); -ProjDataU16 DownsampleProjData(const ProjDataU16& src_proj, const CoordScalar ds_factor); -ProjDataU8 DownsampleProjData(const ProjDataU8& src_proj, const CoordScalar ds_factor); - -ProjDataF32List DownsampleProjData(const ProjDataF32List& src_projs, const CoordScalar ds_factor); -ProjDataU16List DownsampleProjData(const ProjDataU16List& src_projs, const CoordScalar ds_factor); -ProjDataU8List DownsampleProjData(const ProjDataU8List& src_projs, const CoordScalar ds_factor); +ProjDataF32 DownsampleProjData(const ProjDataF32& src_proj, const CoordScalar ds_factor, + const bool force_even_dims = false); +ProjDataU16 DownsampleProjData(const ProjDataU16& src_proj, const CoordScalar ds_factor, + const bool force_even_dims = false); +ProjDataU8 DownsampleProjData(const ProjDataU8& src_proj, const CoordScalar ds_factor, + const bool force_even_dims = false); + +ProjDataF32List DownsampleProjData(const ProjDataF32List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims = false); +ProjDataU16List DownsampleProjData(const ProjDataU16List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims = false); +ProjDataU8List DownsampleProjData(const ProjDataU8List& src_projs, const CoordScalar ds_factor, + const bool force_even_dims = false); template using CamImgPair = std::tuple::Pointer>; diff --git a/lib/itk/xregITKIOUtils.h b/lib/itk/xregITKIOUtils.h index d720f7c..b6388b7 100644 --- a/lib/itk/xregITKIOUtils.h +++ b/lib/itk/xregITKIOUtils.h @@ -121,11 +121,11 @@ typename tImage::Pointer ReadITKImageFromDisk(const std::string& path) return reader->GetOutput(); } -template -typename itk::Image::Pointer -ReadDICOM2DFromDisk(const std::string& path) +template +typename itk::Image::Pointer +ReadDICOMNDFromDisk(const std::string& path) { - using ImgReader = itk::ImageFileReader>; + using ImgReader = itk::ImageFileReader>; typename ImgReader::Pointer img_reader = ImgReader::New(); img_reader->SetFileName(path); @@ -137,6 +137,21 @@ ReadDICOM2DFromDisk(const std::string& path) return img_reader->GetOutput(); } + +template +typename itk::Image::Pointer +ReadDICOM2DFromDisk(const std::string& path) +{ + return ReadDICOMNDFromDisk(path); +} + +template +typename itk::Image::Pointer +ReadDICOM3DFromDisk(const std::string& path) +{ + return ReadDICOMNDFromDisk(path); +} + template void WriteITKLabelMapAsRGB(const itk::Image* img, const std::string& path) { diff --git a/lib/regi/pnp_solvers/xregRANSACPnP.cpp b/lib/regi/pnp_solvers/xregRANSACPnP.cpp index 67b886c..63c3681 100644 --- a/lib/regi/pnp_solvers/xregRANSACPnP.cpp +++ b/lib/regi/pnp_solvers/xregRANSACPnP.cpp @@ -27,67 +27,6 @@ #include "xregAssert.h" #include "xregSampleUtils.h" -std::vector> -xreg::BruteForce3Combos(const size_type num_pts) -{ - std::vector> combos; - - std::vector tmp_arr(3); - - for (size_type i = 0; i < (num_pts - 2); ++i) - { - tmp_arr[0] = i; - - for (size_type j = (i + 1); j < (num_pts - 1); ++j) - { - tmp_arr[1] = j; - - for (size_type k = (j + 1); k < num_pts; ++k) - { - tmp_arr[2] = k; - - combos.push_back(tmp_arr); - } - } - } - - return combos; -} - -std::vector> -xreg::BruteForce4Combos(const size_type num_pts) -{ - xregASSERT(num_pts >= 4); - - std::vector> combos; - - std::vector tmp_arr(4); - - for (size_type i = 0; i < (num_pts - 3); ++i) - { - tmp_arr[0] = i; - - for (size_type j = (i + 1); j < (num_pts - 2); ++j) - { - tmp_arr[1] = j; - - for (size_type k = (j + 1); k < (num_pts - 1); ++k) - { - tmp_arr[2] = k; - - for (size_type l = (k + 1); l < num_pts; ++l) - { - tmp_arr[3] = l; - - combos.push_back(tmp_arr); - } - } - } - } - - return combos; -} - void xreg::RANSACPnP::run_impl() { xregASSERT(pnp_prop_ && pnp_); diff --git a/lib/regi/pnp_solvers/xregRANSACPnP.h b/lib/regi/pnp_solvers/xregRANSACPnP.h index 29fbbc6..26f602f 100644 --- a/lib/regi/pnp_solvers/xregRANSACPnP.h +++ b/lib/regi/pnp_solvers/xregRANSACPnP.h @@ -30,12 +30,6 @@ namespace xreg { -std::vector> -BruteForce3Combos(const size_type num_pts); - -std::vector> -BruteForce4Combos(const size_type num_pts); - class RANSACPnP : public Landmark2D3DRegi { public: diff --git a/lib/spatial/xregFitCircle.cpp b/lib/spatial/xregFitCircle.cpp index 85b9106..0588908 100644 --- a/lib/spatial/xregFitCircle.cpp +++ b/lib/spatial/xregFitCircle.cpp @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020,2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -25,6 +25,7 @@ #include "xregFitCircle.h" #include "xregCMAESInterface.h" +#include "xregSampleUtils.h" std::tuple xreg::FitCircle2D(const Pt2& x, @@ -161,3 +162,103 @@ std::tuple xreg::FitCircle2D(const Pt2List& pts) return std::make_tuple(Pt2(end_x.head(2)), end_x(2)); } + +std::tuple +xreg::FitCircle2DRansac(const Pt2List& pts, const int num_proposals, + const CoordScalar inlier_thresh) +{ + const size_type num_pts = pts.size(); + + std::vector> combos; + + { + std::mt19937 rng_eng; + SeedRNGEngWithRandDev(&rng_eng); + + // for a small enough number of points, brute force the combos. + // also brute force when num_proposals indicates that all combinations should be used + if ((num_proposals <= 0) || (num_pts < 200)) + { + combos = BruteForce3Combos(num_pts); + + // when num_prosals is less than the maximum number of combinations, randomly select + // the desired amount of combos to proposa + if ((num_proposals > 0) && (static_cast(num_proposals) < combos.size())) + { + std::shuffle(combos.begin(), combos.end(), rng_eng); + + combos.resize(num_proposals); + } + } + else + { + // for a large enough number of proposals and for which we do not need all possible + // proposals, uniformly sample some combos + combos = SampleCombos(num_pts, 3, num_proposals, rng_eng); + } + } + + xregASSERT(!combos.empty()); + + size_type cur_best_num_inliers = 0; + + CoordScalar cur_best_mean_error = std::numeric_limits::max(); + + Pt2 cur_best_center; + CoordScalar cur_best_radius; + + Pt2 tmp_center; + CoordScalar tmp_radius; + + Pt2List tmp_inliers; + tmp_inliers.reserve(num_pts); + + for (const auto& cur_combo : combos) + { + std::tie(tmp_center,tmp_radius) = FitCircle2D(pts[cur_combo[0]], pts[cur_combo[1]], + pts[cur_combo[2]]); + + tmp_inliers.clear(); + + for (const auto& p : pts) + { + const CoordScalar p_error = std::abs((p - tmp_center).norm() - tmp_radius); + + if (p_error < inlier_thresh) + { + tmp_inliers.push_back(p); + } + } + + const size_type num_inliers = tmp_inliers.size(); + + if (num_inliers >= cur_best_num_inliers) + { + // this could be the best solution, recompute the circle using all of the inliers + // and recompute the errors + std::tie(tmp_center,tmp_radius) = FitCircle2D(tmp_inliers); + + CoordScalar tmp_error = 0; + + for (const auto& p : tmp_inliers) + { + tmp_error += std::abs((p - tmp_center).norm() - tmp_radius); + } + + tmp_error /= static_cast(num_inliers); + + // keep this solution as the best if it has more inliers than the current best + // OR has a smaller error than the current best + if ((num_inliers > cur_best_num_inliers) || (tmp_error < cur_best_mean_error)) + { + cur_best_num_inliers = num_inliers; + cur_best_mean_error = tmp_error; + cur_best_center = tmp_center; + cur_best_radius = tmp_radius; + } + } + } + + return std::make_tuple(cur_best_center,cur_best_radius); +} + diff --git a/lib/spatial/xregFitCircle.h b/lib/spatial/xregFitCircle.h index 27dbbdf..a305bb4 100644 --- a/lib/spatial/xregFitCircle.h +++ b/lib/spatial/xregFitCircle.h @@ -1,7 +1,7 @@ /* * MIT License * - * Copyright (c) 2020 Robert Grupp + * Copyright (c) 2020,2021 Robert Grupp * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -48,6 +48,17 @@ FitCircle2D(const Pt2& x, /// Requires at least three points. std::tuple FitCircle2D(const Pt2List& pts); +/// \brief Fits a 2D circle to a collection of points using a RANSAC +/// strategy to detect outliers. +/// +/// The FitCircle2D() call with three points is used to generate candidate +/// solutions and create the consensus sets. The call to FitCircle2D() on a +/// collection of points is called using the conensus to produce the final +/// solution. +std::tuple +FitCircle2DRansac(const Pt2List& pts, const int num_proposals = -1, + const CoordScalar inlier_thresh = 0.01f); + } // xreg #endif diff --git a/lib/transforms/xregPerspectiveXform.cpp b/lib/transforms/xregPerspectiveXform.cpp index 952d5d9..f42d67e 100644 --- a/lib/transforms/xregPerspectiveXform.cpp +++ b/lib/transforms/xregPerspectiveXform.cpp @@ -640,7 +640,8 @@ xreg::Pt3 xreg::CalcSourcePositionDelta(const CameraModel& cam1, const CameraMod return delta_cam1_src; } -xreg::CameraModel xreg::DownsampleCameraModel(const CameraModel& src_cam, const CoordScalar ds_factor) +xreg::CameraModel xreg::DownsampleCameraModel(const CameraModel& src_cam, const CoordScalar ds_factor, + const bool force_even_dims) { CameraModel dst_cam; dst_cam.coord_frame_type = src_cam.coord_frame_type; @@ -651,9 +652,24 @@ xreg::CameraModel xreg::DownsampleCameraModel(const CameraModel& src_cam, const intrins(0,2) *= ds_factor; intrins(1,2) *= ds_factor; + long num_ds_rows = std::lround(src_cam.num_det_rows * ds_factor); + long num_ds_cols = std::lround(src_cam.num_det_cols * ds_factor); + + if (force_even_dims) + { + if (num_ds_rows % 2) + { + --num_ds_rows; + } + + if (num_ds_cols % 2) + { + --num_ds_cols; + } + } + dst_cam.setup(intrins, src_cam.extrins.matrix(), - std::lround(src_cam.num_det_rows * ds_factor), - std::lround(src_cam.num_det_cols * ds_factor), + num_ds_rows, num_ds_cols, src_cam.det_row_spacing / ds_factor, src_cam.det_col_spacing / ds_factor); diff --git a/lib/transforms/xregPerspectiveXform.h b/lib/transforms/xregPerspectiveXform.h index d4a9f9f..158aeab 100644 --- a/lib/transforms/xregPerspectiveXform.h +++ b/lib/transforms/xregPerspectiveXform.h @@ -329,7 +329,8 @@ Pt3 CalcSourcePositionDelta(const CameraModel& cam1, const CameraModel& cam2); /// /// A factor of 1 retains the original size, a factor less than 1 downsamples, /// and a factor greater than 1 upsamples. -CameraModel DownsampleCameraModel(const CameraModel& src_cam, const CoordScalar ds_factor); +CameraModel DownsampleCameraModel(const CameraModel& src_cam, const CoordScalar ds_factor, + const bool force_even_dims = false); /// \brief Create a new camera world (extrinsic) frame based on a collection of /// frame transforms from each camera to the new frame.