diff --git a/src/Synergistic/tests/test_cSynergistic.cpp b/src/Synergistic/tests/test_cSynergistic.cpp index dfcd78fa2..95d9940da 100644 --- a/src/Synergistic/tests/test_cSynergistic.cpp +++ b/src/Synergistic/tests/test_cSynergistic.cpp @@ -64,25 +64,22 @@ int main(int argc, char* argv[]) if (argc > 2) mr_recon_h5_filename = argv[2]; - // Test STIR -> Nifti + // Test Nifti -> STIR -> Nifti { std::cout << "\nPerforming STIRImageData to NiftiImageData conversion...\n"; // Load the image as a NiftiImageData3D NiftiImageData3D image_nifti(nifti_filename); - // Read as STIRImageData, convert NiftiImageData3D and save to file - STIRImageData image_stir(nifti_filename); - NiftiImageData3D image_nifti_from_stir(image_stir); - image_nifti_from_stir.write("results/stir_to_nifti.nii",image_nifti.get_original_datatype()); + // Convert to STIR + STIRImageData image_stir_from_nifti(image_nifti); - // Compare the two - if (image_nifti != image_nifti_from_stir) - throw std::runtime_error("Conversion from STIR to Nifti failed"); + // Convert back to NiftiImageData3D + NiftiImageData3D image_nifti_from_stir_from_nifti(image_stir_from_nifti); - // Also save the STIRImageData to file (might be useful visual for comparison) - create_stir_output_file_format("results/stir_output_file_format_nifti.par"); - image_stir.write("results/stir.nii","results/stir_output_file_format_nifti.par"); + // Compare Nifti's + if (image_nifti != image_nifti_from_stir_from_nifti) + throw std::runtime_error("Conversion Nifti->STIR->Nifti failed"); } // Test Gadgetron -> Nifti diff --git a/src/common/GeometricalInfo.cpp b/src/common/GeometricalInfo.cpp index 0d6a2dfcb..db69bff94 100644 --- a/src/common/GeometricalInfo.cpp +++ b/src/common/GeometricalInfo.cpp @@ -20,6 +20,7 @@ limitations under the License. #include "sirf/common/GeometricalInfo.h" #include +#include using namespace sirf; @@ -48,6 +49,87 @@ print_info() const std::cout << "\n"; } +template +typename VoxelisedGeometricalInfo::Coordinate +VoxelisedGeometricalInfo:: +multiply_by_direction_matrix(const DirectionMatrix &matrix, const Coordinate &input) +{ + Coordinate output = {0.f, 0.f, 0.f}; + for (unsigned i=0; i<3; ++i) + for (unsigned j=0; j<3; ++j) + output[i] += matrix[i][j] * input[j]; + return output; +} + +template +typename VoxelisedGeometricalInfo::Index +VoxelisedGeometricalInfo:: +multiply_by_direction_matrix(const DirectionMatrix &matrix, const Index &input) +{ + Coordinate temp = {0.f, 0.f, 0.f}; + for (unsigned i=0; i<3; ++i) + for (unsigned j=0; j<3; ++j) + temp[i] += matrix[i][j] * input[j]; + Index output = { unsigned(temp[0]),unsigned(temp[1]),unsigned(temp[2]) }; + return output; +} + +template +typename VoxelisedGeometricalInfo::DirectionMatrix +VoxelisedGeometricalInfo:: +multiply_direction_matrices(const DirectionMatrix &input1, const DirectionMatrix &input2) +{ + DirectionMatrix output; + for (unsigned i=0; i<3; ++i) + for (unsigned j=0; j<3; ++j) + output[i][j] = 0.f; + for (unsigned i=0; i<3; ++i) + for (unsigned j=0; j<3; ++j) + for (unsigned k=0; k<3; ++k) + output[i][j] += input1[i][k] * input2[k][j]; + + return output; +} + +template +typename VoxelisedGeometricalInfo::DirectionMatrix +VoxelisedGeometricalInfo:: +inverse_direction_matrix(const DirectionMatrix &dm) +{ + DirectionMatrix output; + + // computes the inverse of a matrix m + float det = dm[0][0] * (dm[1][1] * dm[2][2] - dm[2][1] * dm[1][2]) - + dm[0][1] * (dm[1][0] * dm[2][2] - dm[1][2] * dm[2][0]) + + dm[0][2] * (dm[1][0] * dm[2][1] - dm[1][1] * dm[2][0]); + + float invdet = 1.f / det; + + output[0][0] = (dm[1][1] * dm[2][2] - dm[2][1] * dm[1][2]) * invdet; + output[0][1] = (dm[0][2] * dm[2][1] - dm[0][1] * dm[2][2]) * invdet; + output[0][2] = (dm[0][1] * dm[1][2] - dm[0][2] * dm[1][1]) * invdet; + output[1][0] = (dm[1][2] * dm[2][0] - dm[1][0] * dm[2][2]) * invdet; + output[1][1] = (dm[0][0] * dm[2][2] - dm[0][2] * dm[2][0]) * invdet; + output[1][2] = (dm[1][0] * dm[0][2] - dm[0][0] * dm[1][2]) * invdet; + output[2][0] = (dm[1][0] * dm[2][1] - dm[2][0] * dm[1][1]) * invdet; + output[2][1] = (dm[2][0] * dm[0][1] - dm[0][0] * dm[2][1]) * invdet; + output[2][2] = (dm[0][0] * dm[1][1] - dm[1][0] * dm[0][1]) * invdet; + + return output; +} + +template +typename VoxelisedGeometricalInfo::DirectionMatrix +VoxelisedGeometricalInfo:: +absolute_direction_matrix(const DirectionMatrix &dm) +{ + DirectionMatrix output; + for (unsigned i=0; i<4; ++i) + for (unsigned j=0; j<4; ++j) + output[i][j] = fabs(dm[i][j]); + return output; +} + template VoxelisedGeometricalInfo:: VoxelisedGeometricalInfo( diff --git a/src/common/include/sirf/common/GeometricalInfo.h b/src/common/include/sirf/common/GeometricalInfo.h index 21119ef33..20deb8575 100644 --- a/src/common/include/sirf/common/GeometricalInfo.h +++ b/src/common/include/sirf/common/GeometricalInfo.h @@ -91,6 +91,17 @@ class VoxelisedGeometricalInfo : /// Print info virtual void print_info() const; + /// Multiply coordinate by direction matrix + static Coordinate multiply_by_direction_matrix(const DirectionMatrix &matrix, const Coordinate &input); + /// Multiply index by direction matrix + static Index multiply_by_direction_matrix(const DirectionMatrix &matrix, const Index &input); + /// Multiply direction matrix by direction matrix + static DirectionMatrix multiply_direction_matrices(const DirectionMatrix &input1, const DirectionMatrix &input2); + /// Get inverse of direction matrix + static DirectionMatrix inverse_direction_matrix(const DirectionMatrix &dm); + /// Get absolute of direction matrix + static DirectionMatrix absolute_direction_matrix(const DirectionMatrix &dm); + private: Offset _offset; Spacing _spacing; diff --git a/src/xSTIR/cSTIR/stir_data_containers.cpp b/src/xSTIR/cSTIR/stir_data_containers.cpp index d107099a8..26c72603b 100644 --- a/src/xSTIR/cSTIR/stir_data_containers.cpp +++ b/src/xSTIR/cSTIR/stir_data_containers.cpp @@ -260,27 +260,121 @@ const DataContainer& a_y STIRImageData::STIRImageData(const ImageData& id) { - throw std::runtime_error("TODO - create STIRImageData from general SIRFImageData."); - /* The following is incorrect. - Dimensions dim = id.dimensions(); - int nx = dim["x"]; - int ny = dim["y"]; - int nz = 1; - Dimensions::iterator it = dim.begin(); - while (it != dim.end()) { - if (it->first != "x" && it->first != "y") - nz *= it->second; - ++it; + typedef VoxelisedGeometricalInfo3D GeomInfo; + typedef GeomInfo::Size SIRFSize; + typedef GeomInfo::Spacing SIRFSpacing; + typedef GeomInfo::Offset SIRFOffset; + typedef GeomInfo::DirectionMatrix SIRFDM; + typedef GeomInfo::Index SIRFIndex; + + // Get geometrical info + const std::shared_ptr geom_info_sptr = + id.get_geom_info_sptr(); + + // Construct a temporary STIR ImageData and construct a temporary SIRF ImageData based on this + Voxels3DF temp_stir_image( + stir::IndexRange3D(0,1, 0,1, 0,1), + Coord3DF(0.f, 0.f, 0.f), // offset + Coord3DF(1.f, 1.f, 1.f)); // spacing + STIRImageData temp_sirf_image(temp_stir_image); + + // Get the temporary and input direction matrices, multiply them + const SIRFDM &temp_dm = temp_sirf_image.get_geom_info_sptr()->get_direction(); + const SIRFDM &temp_dm_inv = GeomInfo::inverse_direction_matrix(temp_dm); + const SIRFDM &id_dm = geom_info_sptr->get_direction(); + const SIRFDM composed_tm = GeomInfo::multiply_direction_matrices(temp_dm_inv,id_dm); + + std::cout << "\n temp_dm\n"; + for (unsigned i=0;i<3;++i) { + for (unsigned j=0;j<3;++j) { + std::cout << temp_dm[i][j] << " "; + } + std::cout << "\n"; + } + + std::cout << "\n temp_dm_inv\n"; + for (unsigned i=0;i<3;++i) { + for (unsigned j=0;j<3;++j) { + std::cout << temp_dm_inv[i][j] << " "; + } + std::cout << "\n"; + } + + std::cout << "\n id_dm\n"; + for (unsigned i=0;i<3;++i) { + for (unsigned j=0;j<3;++j) { + std::cout << id_dm[i][j] << " "; + } + std::cout << "\n"; + } + + std::cout << "\n composed_tm\n"; + for (unsigned i=0;i<3;++i) { + for (unsigned j=0;j<3;++j) { + std::cout << composed_tm[i][j] << " "; + } + std::cout << "\n"; + } + + const SIRFDM abs_composed_tm = GeomInfo::absolute_direction_matrix(composed_tm); + + std::cout << "\n composed_tm\n"; + for (unsigned i=0;i<3;++i) { + for (unsigned j=0;j<3;++j) { + std::cout << abs_composed_tm[i][j] << " "; + } + std::cout << "\n"; + } + + // Get the input size and spacing, and apply the transformation matrix + const SIRFSize &id_size = geom_info_sptr->get_size(); + const SIRFSpacing &id_spacing = geom_info_sptr->get_spacing(); + const SIRFSize stir_size = GeomInfo::multiply_by_direction_matrix(abs_composed_tm, id_size); + const SIRFSpacing stir_spacing = GeomInfo::multiply_by_direction_matrix(abs_composed_tm, id_spacing); + + Coord3DI size = {int(stir_size[2]),int(stir_size[1]),int(stir_size[0])}; + Coord3DF spacing = {stir_spacing[2], stir_spacing[1], stir_spacing[0]}; + + const Coord3DI stir_min_index = { 0, -(size.y() / 2), -(size.x() / 2) }; + const Coord3DI stir_max_index = stir_min_index + size - 1; + + // From the correct size, get the correct index range + stir::IndexRange<3> index_range(stir_min_index,stir_max_index); + + // Create the actual STIR image + _data = std::make_shared( + index_range, + Coord3DF(0.f, 0.f, 0.f), // offset, will calculate + spacing); + + // Get SIRF offset (as this will be first voxel) + const SIRFOffset &offset = geom_info_sptr->get_offset(); + + const Coord3DF C1 = _data->get_physical_coordinates_for_LPS_coordinates( + _data->get_LPS_coordinates_for_indices(stir_min_index)); + const Coord3DF C2 = _data->get_physical_coordinates_for_LPS_coordinates( + Coord3DF(offset[0],offset[1],offset[2])); + + _data->set_origin(_data->get_origin()-(C1-C2)); + + // Fill the data + GeomInfo::Index idx_id, idx_stir; + for (idx_id[0]=0; idx_id[0]set_up_geom_info();*/ + this->set_up_geom_info(); } void