Implemented safer C++ DataContainer data type handling, fixes #1193

src/Registration/cReg/NiftiBasedRegistration.cpp

@@ -35,7 +35,7 @@ limitations under the License.
 using namespace sirf;
 
 template<class dataType>
-void NiftiBasedRegistration<dataType>::convert_to_NiftiImageData_if_not_already(std::shared_ptr<const NiftiImageData3D<dataType> > &output_sptr, const std::shared_ptr<const ImageData> &input_sptr)
+void NiftiBasedRegistration<dataType>::convert_to_NiftiImageData_if_not_already(std::shared_ptr<const NiftiImageData3D<dataType> > &output_sptr, const std::shared_ptr<const ImageData<dataType> > &input_sptr)
 {
     // Try to dynamic cast from ImageData to (const) NiftiImageData. This will only succeed if original type was NiftiImageData
     output_sptr = std::dynamic_pointer_cast<const NiftiImageData3D<dataType> >(input_sptr);

src/Registration/cReg/NiftiImageData.cpp

@@ -64,27 +64,72 @@ static void check_folder_exists_if_not_create(const std::string &path)
 template<class dataType>
 NiftiImageData<dataType>::NiftiImageData(const NiftiImageData<dataType>& to_copy)
 {
-    *this = dynamic_cast<const ImageData&>(to_copy);
+//std::cout << "in NiftiImageData(const NiftiImageData<dataType>& to_copy)...\n";
+    *this = dynamic_cast<const ImageData<dataType>&>(to_copy);
 }
 
 template<class dataType>
 NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const NiftiImageData<dataType>& to_copy)
 {
-    *this = dynamic_cast<const ImageData&>(to_copy);
+    *this = dynamic_cast<const ImageData<dataType>&>(to_copy);
     return *this;
 }
 
 template<class dataType>
-NiftiImageData<dataType>::NiftiImageData(const ImageData& to_copy)
+NiftiImageData<dataType>::NiftiImageData(const ImageData<float>& to_copy)
 {
+//std::cout << "in NiftiImageData(const ImageData<dataType>& to_copy)...\n";
     *this = to_copy;
 }
 
 template<class dataType>
-NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const ImageData& to_copy)
+NiftiImageData<dataType>::NiftiImageData(const ImageData<complex_float_t>& to_copy)
+{
+//std::cout << "in NiftiImageData(const ImageData<dataType>& to_copy)...\n";
+    *this = to_copy;
+}
+
+template<class dataType>
+NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const ImageData<float>& to_copy)
+{
+//std::cout << "in operator=...\n";
+    // Check for self-assignment
+    if ((void*)this != (void*)&to_copy) {
+//std::cout << "here...\n";
+        // Try to cast to NiftiImageData.
+        const NiftiImageData<dataType> * const nii_ptr = dynamic_cast<const NiftiImageData<dataType> * const >(&to_copy);
+        if (nii_ptr) {
+            // Check the image is copyable
+            if (!nii_ptr->is_initialised())
+                throw std::runtime_error("Trying to copy an uninitialised image.");
+
+            copy_nifti_image(_nifti_image,nii_ptr->_nifti_image);
+            this->_data = static_cast<float*>(_nifti_image->data);
+            this->_original_datatype = nii_ptr->_original_datatype;
+            set_up_geom_info();
+        }
+        else {
+//std::cout << "there...\n";
+            this->_nifti_image = NiftiImageData<float>::create_from_geom_info(*to_copy.get_geom_info_sptr());
+            // Always float
+            this->set_up_data(NIFTI_TYPE_FLOAT32);
+            // Finally, copy the data
+            auto &it_src = to_copy.begin();
+            auto &it_dst = this->begin();
+            for (; it_src != to_copy.end(); ++it_src, ++it_dst)
+                *it_dst = *it_src;
+            set_up_geom_info();
+//std::cout << "ok...\n";
+        }
+    }
+    return *this;
+}
+
+template<class dataType>
+NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const ImageData<complex_float_t>& to_copy)
 {
     // Check for self-assignment
-    if (this != &to_copy) {
+    if ((void*)this != (void*)&to_copy) {
         // Try to cast to NiftiImageData.
         const NiftiImageData<dataType> * const nii_ptr = dynamic_cast<const NiftiImageData<dataType> * const >(&to_copy);
         if (nii_ptr) {
@@ -98,11 +143,16 @@ NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const ImageData& t
             set_up_geom_info();
         }
         else {
+        //std::cout << "here\n" << std::flush;
             this->_nifti_image = NiftiImageData<float>::create_from_geom_info(*to_copy.get_geom_info_sptr());
             // Always float
             this->set_up_data(NIFTI_TYPE_FLOAT32);
             // Finally, copy the data
-            this->copy(to_copy.begin(), this->begin(), this->end());
+            auto &it_src = to_copy.begin();
+            auto &it_dst = this->begin();
+            for (; it_src != to_copy.end(); ++it_src, ++it_dst)
+                *it_dst = (*it_src).complex_float().real();
+            set_up_geom_info();
         }
     }
     return *this;
@@ -213,20 +263,26 @@ std::shared_ptr<nifti_image> NiftiImageData<dataType>::create_from_geom_info(con
 }
 
 template<class dataType>
-void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr<NiftiImageData> &out_sptr, const std::shared_ptr<const ImageData> in_sptr)
+void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr<NiftiImageData>& out_sptr, const std::shared_ptr<const DataContainer> dc_sptr)
 {
-    // Create image from input
+//std::cout << "in construct_NiftiImageData_from_complex_im_real_component...\n";
+    auto in_sptr = std::dynamic_pointer_cast<const ImageData<complex_float_t> >(dc_sptr);
     out_sptr = std::make_shared<NiftiImageData<dataType> >(*in_sptr);
 
     auto &it_in = in_sptr->begin();
     auto &it_out = out_sptr->begin();
     for (; it_in!=in_sptr->end(); ++it_in, ++it_out)
         *it_out = (*it_in).complex_float().real();
+//    out_sptr->set_up_data(NIFTI_TYPE_FLOAT32);
+//    out_sptr->set_up_geom_info();
+//std::cout << "leaving...\n";
 }
 
 template<class dataType>
-void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr<NiftiImageData> &out_sptr, const std::shared_ptr<const ImageData> in_sptr)
+void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr<NiftiImageData>& out_sptr, const std::shared_ptr<const DataContainer> dc_sptr)
 {
+//std::cout << "in construct_NiftiImageData_from_complex_im_imag_component...\n";
+    auto in_sptr = std::dynamic_pointer_cast<const ImageData<complex_float_t> >(dc_sptr);
     if (!in_sptr->is_complex())
         std::cout << "\nNiftiImageData<dataType>::construct_NiftiImageData_from_complex_im. Warning, input image is not complex. Complex component will be empty\n";
 
@@ -237,10 +293,13 @@ void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im_imag_com
     auto &it_out = out_sptr->begin();
     for (; it_in!=in_sptr->end(); ++it_in, ++it_out)
         *it_out = (*it_in).complex_float().imag();
+//    out_sptr->set_up_data(NIFTI_TYPE_FLOAT32);
+//    out_sptr->set_up_geom_info();
+//std::cout << "leaving...\n";
 }
 
 template<class dataType>
-void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im(std::shared_ptr<NiftiImageData> &out_real_sptr, std::shared_ptr<NiftiImageData> &out_imag_sptr, const std::shared_ptr<const ImageData> in_sptr)
+void NiftiImageData<dataType>::construct_NiftiImageData_from_complex_im(std::shared_ptr<NiftiImageData> &out_real_sptr, std::shared_ptr<NiftiImageData> &out_imag_sptr, const std::shared_ptr<const DataContainer> in_sptr)
 {
     construct_NiftiImageData_from_complex_im_real_component(out_real_sptr,in_sptr);
     construct_NiftiImageData_from_complex_im_imag_component(out_imag_sptr,in_sptr);
@@ -394,23 +453,30 @@ void NiftiImageData<dataType>::write(const std::string &filename, const int data
 
     // Create a deep copy in case we need to change datatype
     NiftiImageData copy = *this;
+    //std::cout << "\nwrite ok1\n" << std::flush;
 
     // If user wants to save in a different datatype
     if (datatype != -1)
         copy.change_datatype(datatype);
     else
         copy.change_datatype(_original_datatype);
+    //std::cout << "write ok2\n" << std::flush;
 
     // If ndim==1,but pixdim[2] != 1, set image to 2D (could be 2D image with 1 voxel in 2nd dimension)
     // Similarly, if ndim==2,but pixdim[3] != 1, set image to 3D (could be 3D image with 1 voxel in 3rd dimension)
     for (int i=1; i<=2; ++i)
         if (copy.get_raw_nifti_sptr()->ndim==i && std::abs(copy.get_raw_nifti_sptr()->pixdim[i+1]-1.f) > 1e-4f)
             copy.get_raw_nifti_sptr()->ndim=i+1;
+    //std::cout << "write ok3\n" << std::flush;
+    //std::cout << copy.get_raw_nifti_sptr()->ndim << '\n' << std::flush;
+    //std::cout << copy.get_raw_nifti_sptr()->dim[0] << '\n' << std::flush;
     copy.get_raw_nifti_sptr()->dim[0]=copy.get_raw_nifti_sptr()->ndim;
+    //std::cout << "write ok4\n" << std::flush;
 
     nifti_set_filenames(copy._nifti_image.get(), filename.c_str(), 0, 0);
+    //std::cout << "write ok5\n" << std::flush;
     nifti_image_write(copy._nifti_image.get());
-    std::cout << "done.\n\n";
+    std::cout << "done.\n\n" << std::flush;
 }
 
 template<class dataType>
@@ -1670,10 +1736,11 @@ dataType
 NiftiImageData<dataType>::
 get_inner_product(const NiftiImageData &other) const
 {
-	dataType s;
-	this->dot(other, &s);
-	return s;
-//    return std::inner_product(this->begin(),this->end(),other.begin(),dataType(0));
+    return (dataType)this->dot(other);
+//	dataType s;
+//	this->dot(other, &s);
+//	return s;
+////    return std::inner_product(this->begin(),this->end(),other.begin(),dataType(0));
 }
 
 template<class dataType>
@@ -1738,38 +1805,42 @@ bool NiftiImageData<dataType>::are_equal_to_given_accuracy(const NiftiImageData
 // Pure virtual methods from ImageData
 // ------------------------------------------------------------------------------ //
 template<class dataType>
-void NiftiImageData<dataType>::dot(const DataContainer& a_x, void* ptr) const
+float 
+NiftiImageData<dataType>::dot(const DataContainer& a_x) const
 {
     const NiftiImageData<dataType>& x = dynamic_cast<const NiftiImageData<dataType>&>(a_x);
     ASSERT(_nifti_image->nvox == x._nifti_image->nvox, "dot operands size mismatch");
     double s = 0.0;
     for (unsigned i = 0; i < this->_nifti_image->nvox; ++i)
         s += double(_data[i]) * x._data[i];
-    float* ptr_s = static_cast<float*>(ptr);
-    *ptr_s = (float)s;
+    return (float)s;
+    //float* ptr_s = static_cast<float*>(ptr);
+    //*ptr_s = (float)s;
 }
 
 template<class dataType>
-void NiftiImageData<dataType>::sum(void* ptr) const
+float NiftiImageData<dataType>::sum() const
 {
     double s = 0.0;
     for (unsigned i = 0; i < this->_nifti_image->nvox; ++i)
         s += _data[i];
-    float* ptr_s = static_cast<float*>(ptr);
-    *ptr_s = (float)s;
+    return (float)s;
+    //float* ptr_s = static_cast<float*>(ptr);
+    //*ptr_s = (float)s;
 }
 
 template<class dataType>
-void NiftiImageData<dataType>::max(void* ptr) const
+float NiftiImageData<dataType>::max() const
 {
     float s = 0.0;
     for (unsigned i = 0; i < this->_nifti_image->nvox; ++i) {
         float si = _data[i];
         if (si > s)
             s = si;
     }
-    float* ptr_s = static_cast<float*>(ptr);
-    *ptr_s = s;
+    return s;
+    //float* ptr_s = static_cast<float*>(ptr);
+    //*ptr_s = s;
 }
 
 template<class dataType>
@@ -1919,56 +1990,56 @@ template<class dataType>
 void NiftiImageData<dataType>::multiply
     (const DataContainer& a_x, const DataContainer& a_y)
 {
-    binary_op(a_x, a_y, DataContainer::product<dataType>);
+    binary_op(a_x, a_y, DataContainerTempl<float>::product);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::multiply
 (const DataContainer& a_x, const void* a_y)
 {
-    semibinary_op(a_x, a_y, DataContainer::product<dataType>);
+    semibinary_op(a_x, a_y, DataContainerTempl<float>::product);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::add
 (const DataContainer& a_x, const void* a_y)
 {
-    semibinary_op(a_x, a_y, DataContainer::sum<dataType>);
+    semibinary_op(a_x, a_y, DataContainerTempl<float>::sum);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::divide
     (const DataContainer& a_x, const DataContainer& a_y)
 {
-    binary_op(a_x, a_y, DataContainer::ratio<dataType>);
+    binary_op(a_x, a_y, DataContainerTempl<float>::ratio);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::maximum
 (const DataContainer& a_x, const DataContainer& a_y)
 {
-    binary_op(a_x, a_y, DataContainer::maximum<dataType>);
+    binary_op(a_x, a_y, DataContainerTempl<float>::maximum);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::maximum
 (const DataContainer& a_x, const void* a_y)
 {
-    semibinary_op(a_x, a_y, DataContainer::maximum<dataType>);
+    semibinary_op(a_x, a_y, DataContainerTempl<float>::maximum);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::minimum
 (const DataContainer& a_x, const DataContainer& a_y)
 {
-    binary_op(a_x, a_y, DataContainer::minimum<dataType>);
+    binary_op(a_x, a_y, DataContainerTempl<float>::minimum);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::minimum
 (const DataContainer& a_x, const void* a_y)
 {
-    semibinary_op(a_x, a_y, DataContainer::minimum<dataType>);
+    semibinary_op(a_x, a_y, DataContainerTempl<float>::minimum);
 }
 
 template<class dataType>
@@ -2006,13 +2077,13 @@ void NiftiImageData<dataType>::sqrt(const DataContainer& a_x)
 template<class dataType>
 void NiftiImageData<dataType>::sign(const DataContainer& a_x)
 {
-    unary_op(a_x, DataContainer::sign);
+    unary_op(a_x, DataContainerTempl<float>::sign);
 }
 
 template<class dataType>
 void NiftiImageData<dataType>::abs(const DataContainer& a_x)
 {
-    unary_op(a_x, DataContainer::abs);
+    unary_op(a_x, DataContainerTempl<float>::abs);
 }
 
 template<class dataType>

src/Registration/cReg/NiftiImageData3DTensor.cpp

@@ -181,7 +181,7 @@ template<class dataType>
 void NiftiImageData3DTensor<dataType>::
 tensor_component_maths(
         const int dim,
-        const std::shared_ptr<const ImageData> &scalar_im_sptr,
+        const std::shared_ptr<const ImageData<dataType> > &scalar_im_sptr,
         const typename NiftiImageData<dataType>::MathsType maths_type)
 {
     // Check the dimension to multiply, that dims==5 and nu==3
@@ -218,7 +218,7 @@ tensor_component_maths(
 template<class dataType>
 void NiftiImageData3DTensor<dataType>::
 multiply_tensor_component
-(const int dim, const std::shared_ptr<const ImageData> &scalar_im_sptr)
+(const int dim, const std::shared_ptr<const ImageData<dataType> > &scalar_im_sptr)
 {
     this->tensor_component_maths(dim, scalar_im_sptr, NiftiImageData<dataType>::mul);
 }
@@ -227,7 +227,7 @@ multiply_tensor_component
 template<class dataType>
 void NiftiImageData3DTensor<dataType>::
 add_to_tensor_component
-(const int dim, const std::shared_ptr<const ImageData> &scalar_im_sptr)
+(const int dim, const std::shared_ptr<const ImageData<dataType> > &scalar_im_sptr)
 {
     this->tensor_component_maths(dim, scalar_im_sptr, NiftiImageData<dataType>::ADD);
 }