stir from nifti

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<float> image_nifti(nifti_filename);
 
-            // Read as STIRImageData, convert NiftiImageData3D and save to file
-            STIRImageData image_stir(nifti_filename);
-            NiftiImageData3D<float> 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<float> 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

src/common/GeometricalInfo.cpp

@@ -20,6 +20,7 @@ limitations under the License.
 
 #include "sirf/common/GeometricalInfo.h"
 #include <iostream>
+#include <math.h>
 
 using namespace sirf;
 
@@ -48,6 +49,87 @@ print_info() const
     std::cout << "\n";
 }
 
+template <int num_dimensions>
+typename VoxelisedGeometricalInfo<num_dimensions>::Coordinate
+VoxelisedGeometricalInfo<num_dimensions>::
+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 <int num_dimensions>
+typename VoxelisedGeometricalInfo<num_dimensions>::Index
+VoxelisedGeometricalInfo<num_dimensions>::
+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 <int num_dimensions>
+typename VoxelisedGeometricalInfo<num_dimensions>::DirectionMatrix
+VoxelisedGeometricalInfo<num_dimensions>::
+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 <int num_dimensions>
+typename VoxelisedGeometricalInfo<num_dimensions>::DirectionMatrix
+VoxelisedGeometricalInfo<num_dimensions>::
+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 <int num_dimensions>
+typename VoxelisedGeometricalInfo<num_dimensions>::DirectionMatrix
+VoxelisedGeometricalInfo<num_dimensions>::
+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 <int num_dimensions>
 VoxelisedGeometricalInfo<num_dimensions>::
 VoxelisedGeometricalInfo(

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;

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<const GeomInfo> 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<Voxels3DF>(
+                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]<id_size[0]; ++idx_id[0]) {
+        for (idx_id[1]=0; idx_id[1]<id_size[1]; ++idx_id[1]) {
+            for (idx_id[2]=0; idx_id[2]<id_size[2]; ++idx_id[2]) {
+
+                // Convert input index to stir index
+                idx_stir = GeomInfo::multiply_by_direction_matrix(abs_composed_tm, idx_id);
+
+                // PROBLEM, YOU CAN'T ACCESS IMAGEDATA ELEMENT FROM INDEX.
+                // HOW ABOUT USING FILL AND THEN REORIENTING SOMEHOW?
+                // FUNCTIONALITY MUST ALREADY EXIST IN STIR.
+            }
+        }
     }
-    Voxels3DF voxels(stir::IndexRange3D(0, nz - 1,
-        -(ny / 2), -(ny / 2) + ny - 1, -(nx / 2), -(nx / 2) + nx - 1),
-        Coord3DF(0, 0, 0),
-        Coord3DF(3, 3, 3.375));
-    _data.reset(voxels.clone());
-    copy(id.begin(), begin(), end());
 
     // Set up the geom info
-    this->set_up_geom_info();*/
+    this->set_up_geom_info();
 }
 
 void