diff --git a/src/Registration/cReg/NiftiBasedRegistration.cpp b/src/Registration/cReg/NiftiBasedRegistration.cpp index ff446d9cc..d41d3e252 100644 --- a/src/Registration/cReg/NiftiBasedRegistration.cpp +++ b/src/Registration/cReg/NiftiBasedRegistration.cpp @@ -35,7 +35,7 @@ limitations under the License. using namespace sirf; template -void NiftiBasedRegistration::convert_to_NiftiImageData_if_not_already(std::shared_ptr > &output_sptr, const std::shared_ptr &input_sptr) +void NiftiBasedRegistration::convert_to_NiftiImageData_if_not_already(std::shared_ptr > &output_sptr, const std::shared_ptr > &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 >(input_sptr); diff --git a/src/Registration/cReg/NiftiImageData.cpp b/src/Registration/cReg/NiftiImageData.cpp index 308078984..9557aafd4 100644 --- a/src/Registration/cReg/NiftiImageData.cpp +++ b/src/Registration/cReg/NiftiImageData.cpp @@ -64,27 +64,64 @@ static void check_folder_exists_if_not_create(const std::string &path) template NiftiImageData::NiftiImageData(const NiftiImageData& to_copy) { - *this = dynamic_cast(to_copy); + *this = dynamic_cast&>(to_copy); } template NiftiImageData& NiftiImageData::operator=(const NiftiImageData& to_copy) { - *this = dynamic_cast(to_copy); + *this = dynamic_cast&>(to_copy); return *this; } template -NiftiImageData::NiftiImageData(const ImageData& to_copy) +NiftiImageData::NiftiImageData(const ImageData& to_copy) { *this = to_copy; } template -NiftiImageData& NiftiImageData::operator=(const ImageData& to_copy) +NiftiImageData::NiftiImageData(const ImageData& to_copy) +{ + *this = to_copy; +} + +template +NiftiImageData& NiftiImageData::operator=(const ImageData& to_copy) +{ + // Check for self-assignment + if (static_cast(this) != static_cast(&to_copy)) { + // Try to cast to NiftiImageData. + const NiftiImageData * const nii_ptr = dynamic_cast * 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(_nifti_image->data); + this->_original_datatype = nii_ptr->_original_datatype; + set_up_geom_info(); + } + else { + this->_nifti_image = NiftiImageData::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(); + } + } + return *this; +} + +template +NiftiImageData& NiftiImageData::operator=(const ImageData& to_copy) { // Check for self-assignment - if (this != &to_copy) { + if (static_cast(this) != static_cast(&to_copy)) { // Try to cast to NiftiImageData. const NiftiImageData * const nii_ptr = dynamic_cast * const >(&to_copy); if (nii_ptr) { @@ -102,7 +139,11 @@ NiftiImageData& NiftiImageData::operator=(const ImageData& t // 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,9 +254,9 @@ std::shared_ptr NiftiImageData::create_from_geom_info(con } template -void NiftiImageData::construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr &out_sptr, const std::shared_ptr in_sptr) +void NiftiImageData::construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr& out_sptr, const std::shared_ptr dc_sptr) { - // Create image from input + auto in_sptr = std::dynamic_pointer_cast >(dc_sptr); out_sptr = std::make_shared >(*in_sptr); auto &it_in = in_sptr->begin(); @@ -225,8 +266,9 @@ void NiftiImageData::construct_NiftiImageData_from_complex_im_real_com } template -void NiftiImageData::construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr &out_sptr, const std::shared_ptr in_sptr) +void NiftiImageData::construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr& out_sptr, const std::shared_ptr dc_sptr) { + auto in_sptr = std::dynamic_pointer_cast >(dc_sptr); if (!in_sptr->is_complex()) std::cout << "\nNiftiImageData::construct_NiftiImageData_from_complex_im. Warning, input image is not complex. Complex component will be empty\n"; @@ -240,7 +282,7 @@ void NiftiImageData::construct_NiftiImageData_from_complex_im_imag_com } template -void NiftiImageData::construct_NiftiImageData_from_complex_im(std::shared_ptr &out_real_sptr, std::shared_ptr &out_imag_sptr, const std::shared_ptr in_sptr) +void NiftiImageData::construct_NiftiImageData_from_complex_im(std::shared_ptr &out_real_sptr, std::shared_ptr &out_imag_sptr, const std::shared_ptr 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); @@ -410,7 +452,7 @@ void NiftiImageData::write(const std::string &filename, const int data nifti_set_filenames(copy._nifti_image.get(), filename.c_str(), 0, 0); nifti_image_write(copy._nifti_image.get()); - std::cout << "done.\n\n"; + std::cout << "done.\n\n" << std::flush; } template @@ -1670,10 +1712,7 @@ dataType NiftiImageData:: 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); } template @@ -1738,29 +1777,28 @@ bool NiftiImageData::are_equal_to_given_accuracy(const NiftiImageData // Pure virtual methods from ImageData // ------------------------------------------------------------------------------ // template -void NiftiImageData::dot(const DataContainer& a_x, void* ptr) const +float +NiftiImageData::dot(const DataContainer& a_x) const { const NiftiImageData& x = dynamic_cast&>(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(ptr); - *ptr_s = (float)s; + return (float)s; } template -void NiftiImageData::sum(void* ptr) const +float NiftiImageData::sum() const { double s = 0.0; for (unsigned i = 0; i < this->_nifti_image->nvox; ++i) s += _data[i]; - float* ptr_s = static_cast(ptr); - *ptr_s = (float)s; + return (float)s; } template -void NiftiImageData::max(void* ptr) const +float NiftiImageData::max() const { float s = 0.0; for (unsigned i = 0; i < this->_nifti_image->nvox; ++i) { @@ -1768,10 +1806,60 @@ void NiftiImageData::max(void* ptr) const if (si > s) s = si; } - float* ptr_s = static_cast(ptr); - *ptr_s = s; + return s; } +template +void NiftiImageData::axpby( + dataType a, const DataContainer& a_x, + dataType b, const DataContainer& a_y) +{ + //const float a = *static_cast(ptr_a); + //const float b = *static_cast(ptr_b); + const NiftiImageData& x = dynamic_cast&>(a_x); + const NiftiImageData& y = dynamic_cast&>(a_y); + + // If the result hasn't been initialised, make a clone of one of them + if (!this->is_initialised()) + *this = *x.clone(); + + ASSERT(_nifti_image->nvox == x._nifti_image->nvox, "axpby operands size mismatch"); + ASSERT(_nifti_image->nvox == y._nifti_image->nvox, "axpby operands size mismatch"); + + for (unsigned i=0; i_nifti_image->nvox; ++i) + _data[i] = a * x._data[i] + b * y._data[i]; +} + +template +void NiftiImageData::xapyb( + const DataContainer& a_x, dataType a, + const DataContainer& a_y, dataType b) +{ + NiftiImageData::axpby(a, a_x, b, a_y); +} + +template +void NiftiImageData::xapyb( + const DataContainer& a_x, dataType a, + const DataContainer& a_y, const DataContainer& a_b) +{ + //const float a = *static_cast(ptr_a); + const NiftiImageData& x = dynamic_cast&>(a_x); + const NiftiImageData& y = dynamic_cast&>(a_y); + const NiftiImageData& b = dynamic_cast&>(a_b); + + // If the result hasn't been initialised, make a clone of one of them + if (!this->is_initialised()) + *this = *x.clone(); + + ASSERT(_nifti_image->nvox == x._nifti_image->nvox, "axpby operands size mismatch"); + ASSERT(_nifti_image->nvox == y._nifti_image->nvox, "axpby operands size mismatch"); + ASSERT(_nifti_image->nvox == b._nifti_image->nvox, "axpby operands size mismatch"); + + for (unsigned i = 0; i < this->_nifti_image->nvox; ++i) + _data[i] = a * x._data[i] + b._data[i] * y._data[i]; +} +/* template void NiftiImageData::axpby( const void* ptr_a, const DataContainer& a_x, @@ -1787,9 +1875,9 @@ void NiftiImageData::axpby( *this = *x.clone(); ASSERT(_nifti_image->nvox == x._nifti_image->nvox, "axpby operands size mismatch"); - ASSERT(_nifti_image->nvox == y._nifti_image->nvox, "axpby operands size mismatch"); + ASSERT(_nifti_image->nvox == y._nifti_image->nvox, "axpby operands size mismatch"); - for (unsigned i=0; i_nifti_image->nvox; ++i) + for (unsigned i = 0; i < this->_nifti_image->nvox; ++i) _data[i] = a * x._data[i] + b * y._data[i]; } @@ -1822,7 +1910,7 @@ void NiftiImageData::xapyb( { NiftiImageData::axpby(ptr_a, a_x, ptr_b, a_y); } - +*/ template void NiftiImageData::xapyb( const DataContainer& a_x, const DataContainer& a_a, @@ -1881,11 +1969,12 @@ NiftiImageData::unary_op(const DataContainer& a_x, template void NiftiImageData::semibinary_op(const DataContainer& a_x, - const void* a_y, dataType(*f)(dataType, dataType)) + //const void* a_y, + dataType y, + dataType(*f)(dataType, dataType)) { const NiftiImageData& x = dynamic_cast&>(a_x); - dataType y = *static_cast(a_y); - //dataType y = *(dataType*)a_y; + //dataType y = *static_cast(a_y); // If the result hasn't been initialised, make a clone of one of them if (!this->is_initialised()) @@ -1915,106 +2004,6 @@ void NiftiImageData::binary_op(const DataContainer& a_x, _data[i] = f(x._data[i], y._data[i]); } -template -void NiftiImageData::multiply - (const DataContainer& a_x, const DataContainer& a_y) -{ - binary_op(a_x, a_y, DataContainer::product); -} - -template -void NiftiImageData::multiply -(const DataContainer& a_x, const void* a_y) -{ - semibinary_op(a_x, a_y, DataContainer::product); -} - -template -void NiftiImageData::add -(const DataContainer& a_x, const void* a_y) -{ - semibinary_op(a_x, a_y, DataContainer::sum); -} - -template -void NiftiImageData::divide - (const DataContainer& a_x, const DataContainer& a_y) -{ - binary_op(a_x, a_y, DataContainer::ratio); -} - -template -void NiftiImageData::maximum -(const DataContainer& a_x, const DataContainer& a_y) -{ - binary_op(a_x, a_y, DataContainer::maximum); -} - -template -void NiftiImageData::maximum -(const DataContainer& a_x, const void* a_y) -{ - semibinary_op(a_x, a_y, DataContainer::maximum); -} - -template -void NiftiImageData::minimum -(const DataContainer& a_x, const DataContainer& a_y) -{ - binary_op(a_x, a_y, DataContainer::minimum); -} - -template -void NiftiImageData::minimum -(const DataContainer& a_x, const void* a_y) -{ - semibinary_op(a_x, a_y, DataContainer::minimum); -} - -template -void NiftiImageData::power -(const DataContainer& a_x, const DataContainer& a_y) -{ - binary_op(a_x, a_y, std::pow); -} - -template -void NiftiImageData::power -(const DataContainer& a_x, const void* a_y) -{ - semibinary_op(a_x, a_y, std::pow); -} - -template -void NiftiImageData::exp(const DataContainer& a_x) -{ - unary_op(a_x, std::exp); -} - -template -void NiftiImageData::log(const DataContainer& a_x) -{ - unary_op(a_x, std::log); -} - -template -void NiftiImageData::sqrt(const DataContainer& a_x) -{ - unary_op(a_x, std::sqrt); -} - -template -void NiftiImageData::sign(const DataContainer& a_x) -{ - unary_op(a_x, DataContainer::sign); -} - -template -void NiftiImageData::abs(const DataContainer& a_x) -{ - unary_op(a_x, DataContainer::abs); -} - template void NiftiImageData::set_up_geom_info() { diff --git a/src/Registration/cReg/NiftiImageData3DTensor.cpp b/src/Registration/cReg/NiftiImageData3DTensor.cpp index 532cda4ed..98992ad33 100644 --- a/src/Registration/cReg/NiftiImageData3DTensor.cpp +++ b/src/Registration/cReg/NiftiImageData3DTensor.cpp @@ -181,7 +181,7 @@ template void NiftiImageData3DTensor:: tensor_component_maths( const int dim, - const std::shared_ptr &scalar_im_sptr, + const std::shared_ptr > &scalar_im_sptr, const typename NiftiImageData::MathsType maths_type) { // Check the dimension to multiply, that dims==5 and nu==3 @@ -218,7 +218,7 @@ tensor_component_maths( template void NiftiImageData3DTensor:: multiply_tensor_component -(const int dim, const std::shared_ptr &scalar_im_sptr) +(const int dim, const std::shared_ptr > &scalar_im_sptr) { this->tensor_component_maths(dim, scalar_im_sptr, NiftiImageData::mul); } @@ -227,7 +227,7 @@ multiply_tensor_component template void NiftiImageData3DTensor:: add_to_tensor_component -(const int dim, const std::shared_ptr &scalar_im_sptr) +(const int dim, const std::shared_ptr > &scalar_im_sptr) { this->tensor_component_maths(dim, scalar_im_sptr, NiftiImageData::ADD); } diff --git a/src/Registration/cReg/NiftyResampler.cpp b/src/Registration/cReg/NiftyResampler.cpp index 669323bbf..61ade1dec 100644 --- a/src/Registration/cReg/NiftyResampler.cpp +++ b/src/Registration/cReg/NiftyResampler.cpp @@ -42,18 +42,19 @@ using namespace sirf; using namespace detail; template -static void convert_ImageData_to_ComplexNiftiImageData(ComplexNiftiImageData &output, const std::shared_ptr input_sptr) +static void convert_ImageData_to_ComplexNiftiImageData(ComplexNiftiImageData &output, const std::shared_ptr input_sptr) { // if input is real, only convert first bit if (!input_sptr->is_complex()) { - output.real() = std::make_shared >(*input_sptr); + output.real() = std::make_shared > //(*input_sptr); + (*std::dynamic_pointer_cast >(input_sptr)); output.imag().reset(); } // if input is complex, only set both else { std::shared_ptr > &output_real = output.real(); std::shared_ptr > &output_imag = output.imag(); - NiftiImageData::construct_NiftiImageData_from_complex_im(output_real,output_imag,input_sptr); + NiftiImageData::construct_NiftiImageData_from_complex_im(output_real, output_imag, input_sptr); } } @@ -228,43 +229,49 @@ static void check_images_match( } template -static void set_post_resample_outputs(std::shared_ptr &output_to_return_sptr, std::shared_ptr &output_as_member_sptr, const ComplexNiftiImageData resampled_niftis) +static void set_post_resample_outputs(std::shared_ptr &out_to_return_sptr, std::shared_ptr &output_as_member_sptr, const ComplexNiftiImageData resampled_niftis) { // If output is only real, set that - if (!output_to_return_sptr->is_complex()) + if (!out_to_return_sptr->is_complex()) { + auto output_to_return_sptr = std::dynamic_pointer_cast >(out_to_return_sptr); output_to_return_sptr->fill(*resampled_niftis.real()); + output_as_member_sptr = output_to_return_sptr; + } // Else, set the complex bit else { - NumberType::Type output_num_type = (*output_to_return_sptr->begin()).get_typeID(); - if (output_num_type != NumberType::CXFLOAT) + auto output_to_return_sptr = std::dynamic_pointer_cast >(out_to_return_sptr); + //NumberType::Type output_num_type = (*output_to_return_sptr->begin()).get_typeID(); + //if (output_num_type != NumberType::CXFLOAT) + if (!output_to_return_sptr.get()) throw std::runtime_error("NiftyResampler: Only complex type currently supported is complex float"); - ImageData::Iterator &it_out = output_to_return_sptr->begin(); + typename ImageData ::Iterator &it_out = output_to_return_sptr->begin(); auto &it_real = resampled_niftis.real()->begin(); auto &it_imag = resampled_niftis.imag()->begin(); for (; it_out!=output_to_return_sptr->end(); ++it_real, ++it_imag, ++it_out) { complex_float_t cmplx_flt(*it_real,*it_imag); - *it_out = NumRef((void *)&cmplx_flt, output_num_type); + *it_out = NumRef((void *)&cmplx_flt, NumberType::CXFLOAT); } + output_as_member_sptr = output_to_return_sptr; } // Copy the output so that backwards compatibility of get_output() is preserved. - output_as_member_sptr = output_to_return_sptr; + //output_as_member_sptr = output_to_return_sptr; } template -std::shared_ptr NiftyResampler::forward(const std::shared_ptr input_sptr) +std::shared_ptr NiftyResampler::forward(const std::shared_ptr input_sptr) { // Call the set up set_up_forward(); - std::shared_ptr output_sptr = this->_reference_image_sptr->clone(); + std::shared_ptr output_sptr = this->_reference_image_sptr->clone(); forward(output_sptr, input_sptr); return output_sptr; } template -void NiftyResampler::forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) +void NiftyResampler::forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) { // Call the set up set_up_forward(); @@ -295,19 +302,19 @@ void NiftyResampler::forward(std::shared_ptr output_sptr, c } template -std::shared_ptr NiftyResampler::adjoint(const std::shared_ptr input_sptr) +std::shared_ptr NiftyResampler::adjoint(const std::shared_ptr input_sptr) { // Call the set up set_up_adjoint(); - std::shared_ptr output_sptr = this->_floating_image_sptr->clone(); + std::shared_ptr output_sptr = this->_floating_image_sptr->clone(); adjoint(output_sptr, input_sptr); return output_sptr; } template -void NiftyResampler::adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) +void NiftyResampler::adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) { // Call the set up set_up_adjoint(); diff --git a/src/Registration/cReg/Registration.cpp b/src/Registration/cReg/Registration.cpp index 4e5212690..e053b4050 100644 --- a/src/Registration/cReg/Registration.cpp +++ b/src/Registration/cReg/Registration.cpp @@ -34,14 +34,14 @@ limitations under the License. using namespace sirf; template -void Registration::set_reference_image(const std::shared_ptr reference_image_sptr) +void Registration::set_reference_image(const std::shared_ptr > reference_image_sptr) { _reference_image_sptr = reference_image_sptr; _reference_image_filename = ""; } template -void Registration::set_floating_image(const std::shared_ptr floating_image_sptr) +void Registration::set_floating_image(const std::shared_ptr > floating_image_sptr) { if (!_floating_image_filenames.empty()) { std::cout << "\nClearing floating images set via filename"; @@ -53,7 +53,7 @@ void Registration::set_floating_image(const std::shared_ptr -void Registration::add_floating_image(const std::shared_ptr floating_image_sptr) +void Registration::add_floating_image(const std::shared_ptr > floating_image_sptr) { if (!_floating_image_filenames.empty()) { std::cout << "\nClearing floating images set via filename"; diff --git a/src/Registration/cReg/Resample.cpp b/src/Registration/cReg/Resample.cpp index db50ba0d9..2b965e4dd 100644 --- a/src/Registration/cReg/Resample.cpp +++ b/src/Registration/cReg/Resample.cpp @@ -33,7 +33,7 @@ using namespace sirf; /// Set reference image template -void Resampler::set_reference_image(const std::shared_ptr reference_image_sptr) +void Resampler::set_reference_image(const std::shared_ptr reference_image_sptr) { _reference_image_sptr = reference_image_sptr; _need_to_set_up = true; @@ -43,7 +43,7 @@ void Resampler::set_reference_image(const std::shared_ptr -void Resampler::set_floating_image(const std::shared_ptr floating_image_sptr) +void Resampler::set_floating_image(const std::shared_ptr floating_image_sptr) { _floating_image_sptr = floating_image_sptr; _need_to_set_up = true; @@ -82,13 +82,13 @@ void Resampler::check_parameters() } template -std::shared_ptr Resampler::backward(const std::shared_ptr input_sptr) +std::shared_ptr Resampler::backward(const std::shared_ptr input_sptr) { return adjoint(input_sptr); } template -void Resampler::backward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) +void Resampler::backward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) { adjoint(output_sptr, input_sptr); } diff --git a/src/Registration/cReg/cReg.cpp b/src/Registration/cReg/cReg.cpp index c581671ee..d3c397e39 100644 --- a/src/Registration/cReg/cReg.cpp +++ b/src/Registration/cReg/cReg.cpp @@ -466,12 +466,12 @@ extern "C" void* cReg_NiftiImageData_from_SIRFImageData(void* ptr, const int is_3D) { try { - ImageData& sirf_im = objectFromHandle(ptr); + ImageData& sirf_im = objectFromHandle >(ptr); if (is_3D) - return newObjectHandle(std::make_shared>(sirf_im)); + return newObjectHandle(std::make_shared >(sirf_im)); else - return newObjectHandle(std::make_shared>(sirf_im)); + return newObjectHandle(std::make_shared >(sirf_im)); } CATCH; } @@ -480,8 +480,8 @@ extern "C" void* cReg_NiftiImageData_from_complex_ImageData_real_component(void* in_ptr) { try { - std::shared_ptr in_sptr; - getObjectSptrFromHandle(in_ptr, in_sptr); + std::shared_ptr in_sptr; + getObjectSptrFromHandle(in_ptr, in_sptr); std::shared_ptr > out_sptr; NiftiImageData::construct_NiftiImageData_from_complex_im_real_component(out_sptr, in_sptr); return newObjectHandle(out_sptr); @@ -493,8 +493,8 @@ extern "C" void* cReg_NiftiImageData_from_complex_ImageData_imag_component(void* in_ptr) { try { - std::shared_ptr in_sptr; - getObjectSptrFromHandle(in_ptr, in_sptr); + std::shared_ptr in_sptr; + getObjectSptrFromHandle(in_ptr, in_sptr); std::shared_ptr > out_sptr; NiftiImageData::construct_NiftiImageData_from_complex_im_imag_component(out_sptr, in_sptr); return newObjectHandle(out_sptr); @@ -554,9 +554,9 @@ extern "C" void* cReg_NiftiImageData3DTensor_construct_from_3_components(const char* obj, const void *x_ptr, const void *y_ptr, const void *z_ptr) { try { - ImageData& x = objectFromHandle(x_ptr); - ImageData& y = objectFromHandle(y_ptr); - ImageData& z = objectFromHandle(z_ptr); + ImageData& x = objectFromHandle >(x_ptr); + ImageData& y = objectFromHandle >(y_ptr); + ImageData& z = objectFromHandle >(z_ptr); std::shared_ptr > sptr; if (strcmp(obj,"NiftiImageData3DTensor") == 0) @@ -678,8 +678,8 @@ void* cReg_Registration_add_floating(const void* ptr, const void* im_ptr) { try { Registration& reg = objectFromHandle >(ptr); - std::shared_ptr im_sptr; - getObjectSptrFromHandle(im_ptr, im_sptr); + std::shared_ptr > im_sptr; + getObjectSptrFromHandle >(im_ptr, im_sptr); reg.add_floating_image(im_sptr); return new DataHandle; } @@ -853,10 +853,10 @@ void* cReg_NiftyResampler_forward(const void* output_ptr, const void * const inp getObjectSptrFromHandle >(resampler_ptr, resampler_sptr); // Get input and output images - std::shared_ptr input_sptr; - getObjectSptrFromHandle(input_ptr, input_sptr); - std::shared_ptr output_sptr; - getObjectSptrFromHandle(output_ptr, output_sptr); + std::shared_ptr > input_sptr; + getObjectSptrFromHandle >(input_ptr, input_sptr); + std::shared_ptr > output_sptr; + getObjectSptrFromHandle >(output_ptr, output_sptr); // Forward transformation resampler_sptr->forward(output_sptr,input_sptr); @@ -874,10 +874,10 @@ void* cReg_NiftyResampler_adjoint(const void* output_ptr, const void * const inp getObjectSptrFromHandle >(resampler_ptr, resampler_sptr); // Get input and output images - std::shared_ptr input_sptr; - getObjectSptrFromHandle(input_ptr, input_sptr); - std::shared_ptr output_sptr; - getObjectSptrFromHandle(output_ptr, output_sptr); + std::shared_ptr > input_sptr; + getObjectSptrFromHandle >(input_ptr, input_sptr); + std::shared_ptr > output_sptr; + getObjectSptrFromHandle >(output_ptr, output_sptr); // Forward transformation resampler_sptr->adjoint(output_sptr,input_sptr); diff --git a/src/Registration/cReg/cReg_p.cpp b/src/Registration/cReg/cReg_p.cpp index 26199559d..5454f3961 100644 --- a/src/Registration/cReg/cReg_p.cpp +++ b/src/Registration/cReg/cReg_p.cpp @@ -115,15 +115,15 @@ sirf::cReg_NiftiImageDataParameter(const DataHandle* handle, const char* name) void* sirf::cReg_setRegistrationParameter(void* hp, const char* name, const void* hv) { - std::shared_ptr im_sptr; + std::shared_ptr > im_sptr; Registration& s = objectFromHandle >(hp); if (strcmp(name, "reference_image") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_reference_image(im_sptr); } else if (strcmp(name, "floating_image") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_floating_image(im_sptr); } else @@ -137,17 +137,17 @@ sirf::cReg_setRegistrationParameter(void* hp, const char* name, const void* hv) void* sirf::cReg_setNiftyRegistrationParameter(void* hp, const char* name, const void* hv) { - std::shared_ptr im_sptr; + std::shared_ptr > im_sptr; NiftyRegistration& s = objectFromHandle >(hp); if (strcmp(name, "parameter_file") == 0) s.set_parameter_file(charDataFromHandle(hv)); else if (strcmp(name, "reference_mask") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_reference_mask(im_sptr); } else if (strcmp(name, "floating_mask") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_floating_mask(im_sptr); } else @@ -202,15 +202,15 @@ sirf::cReg_setSPMRegistrationParameter(void* hp, const char* name, const void* h void* sirf::cReg_setNiftyResamplerParameter(void* hp, const char* name, const void* hv) { - std::shared_ptr im_sptr; + std::shared_ptr > im_sptr; NiftyResampler& s = objectFromHandle >(hp); if (strcmp(name, "reference_image") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_reference_image(im_sptr); } else if (strcmp(name, "floating_image") == 0) { - getObjectSptrFromHandle(hv, im_sptr); + getObjectSptrFromHandle >(hv, im_sptr); s.set_floating_image(im_sptr); } else if (strcmp(name, "interpolation_type") == 0) diff --git a/src/Registration/cReg/include/sirf/Reg/NiftiBasedRegistration.h b/src/Registration/cReg/include/sirf/Reg/NiftiBasedRegistration.h index 3d8acc9fd..83f08b55e 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftiBasedRegistration.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftiBasedRegistration.h @@ -61,7 +61,7 @@ class NiftiBasedRegistration : public Registration virtual const std::shared_ptr > get_displacement_field_inverse_sptr(const unsigned idx = 0) const; /// Convert an ImageData to NiftiImageData. Try to dynamic pointer cast, else create new image. - static void convert_to_NiftiImageData_if_not_already(std::shared_ptr > &output_sptr, const std::shared_ptr &input_sptr); + static void convert_to_NiftiImageData_if_not_already(std::shared_ptr > &output_sptr, const std::shared_ptr > &input_sptr); protected: diff --git a/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h b/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h index 09c0f7eae..dbb704400 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h @@ -41,621 +41,621 @@ limitations under the License. namespace sirf { -/*! -\file -\ingroup Registration -\brief Base class for SIRF's nifti image data. + /*! + \file + \ingroup Registration + \brief Base class for SIRF's nifti image data. -This is a wrapper around the basic nifti_image (stored as a sptr), -with extra functionality. + This is a wrapper around the basic nifti_image (stored as a sptr), + with extra functionality. -This is the general form, and any number of dimensions are allowed. -This is contrary to the derived classes, which have specific requirements -(e.g., NiftiImageData3DDeformation requires 5 dimensions, of which the 4th (time) == 1). + This is the general form, and any number of dimensions are allowed. + This is contrary to the derived classes, which have specific requirements + (e.g., NiftiImageData3DDeformation requires 5 dimensions, of which the 4th (time) == 1). -To move between different images types (e.g., STIRImageData and MRImageData), -we need to know the patient's position relative to the scanner. For this, only -qform_code == 1 will suffice. + To move between different images types (e.g., STIRImageData and MRImageData), + we need to know the patient's position relative to the scanner. For this, only + qform_code == 1 will suffice. -qform/sform -The qform code describes "scanner-anatomical" coordinates, whereas the -sform code describes the location of the voxels in some standard space. + qform/sform + The qform code describes "scanner-anatomical" coordinates, whereas the + sform code describes the location of the voxels in some standard space. -For qform > 0, the origin of coordinates would generally be whatever -the scanner origin is; for example, in MRI, (0,0,0) is the center -of the gradient coil. + For qform > 0, the origin of coordinates would generally be whatever + the scanner origin is; for example, in MRI, (0,0,0) is the center + of the gradient coil. -For sform > 0, the origin of coordinates would depend on the value -of sform_code; for example, for the Talairach coordinate system, -(0,0,0) corresponds to the Anterior Commissure. + For sform > 0, the origin of coordinates would depend on the value + of sform_code; for example, for the Talairach coordinate system, + (0,0,0) corresponds to the Anterior Commissure. -\author Richard Brown -\author SyneRBI -*/ + \author Richard Brown + \author SyneRBI + */ -template -class NiftiImageData : public ImageData -{ -public: + //template + //class ImageData; - typedef ImageData::Iterator BaseIter; - typedef ImageData::Iterator_const BaseIter_const; - class Iterator : public BaseIter { - public: - Iterator(dataType *iter) : _iter(iter) - {} - Iterator& operator=(const Iterator& iter) - { - _iter = iter._iter; - _ref.copy(iter._ref); - return *this; - } - virtual Iterator& operator++() - { - ++_iter; - return *this; - } - virtual bool operator==(const BaseIter& an_iter) const - { - const Iterator& iter = (const Iterator&)an_iter; - return _iter == iter._iter; - } - virtual bool operator!=(const BaseIter& an_iter) const - { - const Iterator& iter = (const Iterator&)an_iter; - return _iter != iter._iter; - } - virtual FloatRef& operator*() - { - dataType& v = *_iter; - _ref.set_ptr(&v); - return _ref; - } - private: - dataType *_iter; - FloatRef _ref; - }; - class Iterator_const : public BaseIter_const { + template //, typename datatype = dataType> + class NiftiImageData : public ImageData + { public: - Iterator_const(const dataType *iter) : _iter(iter) - {} - Iterator_const& operator=(const Iterator_const& iter) - { - _iter = iter._iter; - _ref.copy(iter._ref); - return *this; - } - virtual Iterator_const& operator++() - { - ++_iter; - return *this; - } - virtual bool operator==(const BaseIter_const& an_iter) const - { - const Iterator_const& iter = (const Iterator_const&)an_iter; - return _iter == iter._iter; - } - virtual bool operator!=(const BaseIter_const& an_iter) const - { - const Iterator_const& iter = (const Iterator_const&)an_iter; - return _iter != iter._iter; - } - virtual const FloatRef& operator*() const + + typedef typename ImageData::Iterator BaseIter; + typedef typename ImageData::Iterator_const BaseIter_const; + class Iterator : public BaseIter { + public: + Iterator(dataType* iter) : _iter(iter) + {} + Iterator& operator=(const Iterator& iter) + { + _iter = iter._iter; + _ref.copy(iter._ref); + return *this; + } + virtual Iterator& operator++() + { + ++_iter; + return *this; + } + virtual bool operator==(const BaseIter& an_iter) const + { + const Iterator& iter = (const Iterator&)an_iter; + return _iter == iter._iter; + } + virtual bool operator!=(const BaseIter& an_iter) const + { + const Iterator& iter = (const Iterator&)an_iter; + return _iter != iter._iter; + } + virtual FloatRef& operator*() + { + dataType& v = *_iter; + _ref.set_ptr(&v); + return _ref; + } + private: + dataType* _iter; + FloatRef _ref; + }; + class Iterator_const : public BaseIter_const { + public: + Iterator_const(const dataType* iter) : _iter(iter) + {} + Iterator_const& operator=(const Iterator_const& iter) + { + _iter = iter._iter; + _ref.copy(iter._ref); + return *this; + } + virtual Iterator_const& operator++() + { + ++_iter; + return *this; + } + virtual bool operator==(const BaseIter_const& an_iter) const + { + const Iterator_const& iter = (const Iterator_const&)an_iter; + return _iter == iter._iter; + } + virtual bool operator!=(const BaseIter_const& an_iter) const + { + const Iterator_const& iter = (const Iterator_const&)an_iter; + return _iter != iter._iter; + } + virtual const FloatRef& operator*() const + { + const dataType& v = *_iter; + _ref.set_ptr((void*)&v); + return _ref; + } + private: + const dataType* _iter; + mutable FloatRef _ref; + }; + + virtual std::string data_type() const { - const dataType& v = *_iter; - _ref.set_ptr((void*)&v); - return _ref; + return std::string("float"); } - private: - const dataType *_iter; - mutable FloatRef _ref; - }; - /// Constructor - NiftiImageData() {} + /// Constructor + NiftiImageData() {} - /// Destructor - virtual ~NiftiImageData() {} + /// Destructor + virtual ~NiftiImageData() {} - /// Copy constructor - NiftiImageData(const NiftiImageData& to_copy); + /// Copy constructor + NiftiImageData(const NiftiImageData& to_copy); - /// Assignment - NiftiImageData& operator=(const NiftiImageData& to_copy); + /// Assignment + NiftiImageData& operator=(const NiftiImageData& to_copy); - /// Copy constructor - NiftiImageData(const ImageData& to_copy); + /// Copy constructor + NiftiImageData(const DataContainer& to_copy) + { + //std::cout << "in NiftiImageData(const DataContainer& to_copy...)\n"; + const ImageData& im_data = dynamic_cast&>(to_copy); + *this = im_data; + } - /// Assignment - NiftiImageData& operator=(const ImageData& to_copy); + /// Copy constructor + NiftiImageData(const ImageData& to_copy); + NiftiImageData(const ImageData& to_copy); - /// Filename constructor - NiftiImageData(const std::string &filename); + /// Assignment + NiftiImageData& operator=(const ImageData& to_copy); + NiftiImageData& operator=(const ImageData& to_copy); - /// Nifti constructor - NiftiImageData(const nifti_image &image_nifti); + /// Filename constructor + NiftiImageData(const std::string& filename); - /// Construct from array - template - NiftiImageData(const inputType * const data, const VoxelisedGeometricalInfo3D &geom, const bool is_tensor = false) - { - this->_nifti_image = create_from_geom_info(geom, is_tensor); - - // Set the datatype - if (typeid(inputType) == typeid(bool)) this->set_up_data(DT_BINARY); - else if (typeid(inputType) == typeid(signed char)) this->set_up_data(DT_INT8); - else if (typeid(inputType) == typeid(signed short)) this->set_up_data(DT_INT16); - else if (typeid(inputType) == typeid(signed int)) this->set_up_data(DT_INT32); - else if (typeid(inputType) == typeid(float)) this->set_up_data(DT_FLOAT32); - else if (typeid(inputType) == typeid(double)) this->set_up_data(DT_FLOAT64); - else if (typeid(inputType) == typeid(unsigned char)) this->set_up_data(DT_UINT8); - else if (typeid(inputType) == typeid(unsigned short)) this->set_up_data(DT_UINT16); - else if (typeid(inputType) == typeid(unsigned int)) this->set_up_data(DT_UINT32); - else if (typeid(inputType) == typeid(signed long long)) this->set_up_data(DT_INT64); - else if (typeid(inputType) == typeid(unsigned long long)) this->set_up_data(DT_UINT64); - else if (typeid(inputType) == typeid(long double)) this->set_up_data(DT_FLOAT128); - else { - std::stringstream ss; - ss << "NiftiImageData constructor from raw array: "; - ss << typeid (inputType).name(); - ss << " (bytes per voxel: "; - ss << sizeof(inputType) << ")."; - throw std::runtime_error(ss.str()); + /// Nifti constructor + NiftiImageData(const nifti_image& image_nifti); + + /// Construct from array + template + NiftiImageData(const inputType* const data, const VoxelisedGeometricalInfo3D& geom, const bool is_tensor = false) + { + this->_nifti_image = create_from_geom_info(geom, is_tensor); + + // Set the datatype + if (typeid(inputType) == typeid(bool)) this->set_up_data(DT_BINARY); + else if (typeid(inputType) == typeid(signed char)) this->set_up_data(DT_INT8); + else if (typeid(inputType) == typeid(signed short)) this->set_up_data(DT_INT16); + else if (typeid(inputType) == typeid(signed int)) this->set_up_data(DT_INT32); + else if (typeid(inputType) == typeid(float)) this->set_up_data(DT_FLOAT32); + else if (typeid(inputType) == typeid(double)) this->set_up_data(DT_FLOAT64); + else if (typeid(inputType) == typeid(unsigned char)) this->set_up_data(DT_UINT8); + else if (typeid(inputType) == typeid(unsigned short)) this->set_up_data(DT_UINT16); + else if (typeid(inputType) == typeid(unsigned int)) this->set_up_data(DT_UINT32); + else if (typeid(inputType) == typeid(signed long long)) this->set_up_data(DT_INT64); + else if (typeid(inputType) == typeid(unsigned long long)) this->set_up_data(DT_UINT64); + else if (typeid(inputType) == typeid(long double)) this->set_up_data(DT_FLOAT128); + else { + std::stringstream ss; + ss << "NiftiImageData constructor from raw array: "; + ss << typeid (inputType).name(); + ss << " (bytes per voxel: "; + ss << sizeof(inputType) << ")."; + throw std::runtime_error(ss.str()); + } + + for (unsigned i = 0; i < _nifti_image->nvox; ++i) + this->_data[i] = dataType(data[i]); } - for (unsigned i=0; i<_nifti_image->nvox; ++i) - this->_data[i] = dataType(data[i]); - } + /// Create NiftiImageData from geometrical info + static std::shared_ptr create_from_geom_info(const VoxelisedGeometricalInfo3D& geom, const bool is_tensor = false, const NREG_TRANS_TYPE tensor_type = NREG_TRANS_TYPE::DEF_FIELD); + + /// Construct NiftiImageData from the real component of a complex SIRF ImageData + static void construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr& out_sptr, const std::shared_ptr in_sptr); - /// Create NiftiImageData from geometrical info - static std::shared_ptr create_from_geom_info(const VoxelisedGeometricalInfo3D &geom, const bool is_tensor=false, const NREG_TRANS_TYPE tensor_type=NREG_TRANS_TYPE::DEF_FIELD); + /// Construct NiftiImageData from the imaginary component of a complex SIRF ImageData + static void construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr& out_sptr, const std::shared_ptr in_sptr); - /// Construct NiftiImageData from the real component of a complex SIRF ImageData - static void construct_NiftiImageData_from_complex_im_real_component(std::shared_ptr &out_sptr, const std::shared_ptr in_sptr); + /// Construct two NiftiImageData from a complex SIRF ImageData + static void construct_NiftiImageData_from_complex_im(std::shared_ptr& out_real_sptr, std::shared_ptr& out_imag_sptr, const std::shared_ptr in_sptr); - /// Construct NiftiImageData from the imaginary component of a complex SIRF ImageData - static void construct_NiftiImageData_from_complex_im_imag_component(std::shared_ptr &out_sptr, const std::shared_ptr in_sptr); + /// Equality operator + bool operator==(const NiftiImageData& other) const; - /// Construct two NiftiImageData from a complex SIRF ImageData - static void construct_NiftiImageData_from_complex_im(std::shared_ptr &out_real_sptr, std::shared_ptr &out_imag_sptr, const std::shared_ptr in_sptr); + /// Equality operator + bool operator!=(const NiftiImageData& other) const; - /// Equality operator - bool operator==(const NiftiImageData &other) const; + /// Addition operator + NiftiImageData& operator+=(const NiftiImageData& rhs); - /// Equality operator - bool operator!=(const NiftiImageData &other) const; + /// Addition operator + friend NiftiImageData operator+(NiftiImageData lhs, const NiftiImageData& rhs) + { + lhs += rhs; + return lhs; + } - /// Addition operator - NiftiImageData& operator+=(const NiftiImageData &rhs); + /// Subtraction operator + NiftiImageData& operator-=(const NiftiImageData& rhs); - /// Addition operator - friend NiftiImageData operator+(NiftiImageData lhs, const NiftiImageData& rhs) - { - lhs += rhs; - return lhs; - } + /// Subtraction operator + friend NiftiImageData operator-(NiftiImageData lhs, const NiftiImageData& rhs) + { + lhs -= rhs; + return lhs; + } - /// Subtraction operator - NiftiImageData& operator-=(const NiftiImageData &rhs); + /// Addition operator + NiftiImageData& operator+=(const float); - /// Subtraction operator - friend NiftiImageData operator-(NiftiImageData lhs, const NiftiImageData& rhs) - { - lhs -= rhs; - return lhs; - } + /// Addition operator + friend NiftiImageData operator+(NiftiImageData lhs, const float val) + { + lhs += val; + return lhs; + } - /// Addition operator - NiftiImageData& operator+=(const float); + /// Subtraction operator + NiftiImageData& operator-=(const float); - /// Addition operator - friend NiftiImageData operator+(NiftiImageData lhs, const float val) - { - lhs += val; - return lhs; - } + /// Subtraction operator + friend NiftiImageData operator-(NiftiImageData lhs, const float val) + { + lhs -= val; + return lhs; + } - /// Subtraction operator - NiftiImageData& operator-=(const float); + /// Multiplication operator + NiftiImageData& operator*=(const float); + /// Multiplication operator + NiftiImageData& operator*=(const NiftiImageData& rhs); - /// Subtraction operator - friend NiftiImageData operator-(NiftiImageData lhs, const float val) - { - lhs -= val; - return lhs; - } + /// Multiplication operator + friend NiftiImageData operator*(NiftiImageData lhs, const float val) + { + lhs *= val; + return lhs; + } - /// Multiplication operator - NiftiImageData& operator*=(const float); - /// Multiplication operator - NiftiImageData& operator*=(const NiftiImageData &rhs); + /// Multiplication operator + friend NiftiImageData operator*(NiftiImageData lhs, const NiftiImageData& rhs) + { + lhs *= rhs; + return lhs; + } - /// Multiplication operator - friend NiftiImageData operator*(NiftiImageData lhs, const float val) - { - lhs *= val; - return lhs; - } + /// Division operator + NiftiImageData& operator/=(const float); + // /// Division operator + NiftiImageData& operator/=(const NiftiImageData& rhs); - /// Multiplication operator - friend NiftiImageData operator*(NiftiImageData lhs, const NiftiImageData& rhs) - { - lhs *= rhs; - return lhs; - } - - /// Division operator - NiftiImageData& operator/=(const float); - // /// Division operator - NiftiImageData& operator/=(const NiftiImageData &rhs); - - // /// Division operator - friend NiftiImageData operator/(NiftiImageData lhs, const float val) - { - lhs /= val; - return lhs; - } + // /// Division operator + friend NiftiImageData operator/(NiftiImageData lhs, const float val) + { + lhs /= val; + return lhs; + } - /// Division operator - friend NiftiImageData operator/(NiftiImageData lhs, const NiftiImageData& rhs) - { - lhs /= rhs; - return lhs; - } + /// Division operator + friend NiftiImageData operator/(NiftiImageData lhs, const NiftiImageData& rhs) + { + lhs /= rhs; + return lhs; + } - /// Access data element via 1D index (const) - float operator()(const int index) const; + /// Access data element via 1D index (const) + float operator()(const int index) const; - /// Access data element via 1D index - float &operator()(const int index); + /// Access data element via 1D index + float& operator()(const int index); - /// Access data element via 7D index (const) - float operator()(const int index[7]) const; + /// Access data element via 7D index (const) + float operator()(const int index[7]) const; - /// Access data element via 7D index - float &operator()(const int index[7]); + /// Access data element via 7D index + float& operator()(const int index[7]); - /// Access data element via 7D index (const) - float operator()(const int x, const int y, const int z, const int t=0, const int u=0, const int v=0, const int w=0) const; + /// Access data element via 7D index (const) + float operator()(const int x, const int y, const int z, const int t = 0, const int u = 0, const int v = 0, const int w = 0) const; - /// Access data element via 7D index - float &operator()(const int x, const int y, const int z, const int t=0, const int u=0, const int v=0, const int w=0); + /// Access data element via 7D index + float& operator()(const int x, const int y, const int z, const int t = 0, const int u = 0, const int v = 0, const int w = 0); - /// Is the image initialised? (Should be unless default constructor was used.) - bool is_initialised() const { return (_nifti_image && _data && _nifti_image->datatype == DT_FLOAT32 ? true : false); } + /// Is the image initialised? (Should be unless default constructor was used.) + bool is_initialised() const { return (_nifti_image && _data && _nifti_image->datatype == DT_FLOAT32 ? true : false); } - /// Get image as nifti as const - std::shared_ptr get_raw_nifti_sptr() const; + /// Get image as nifti as const + std::shared_ptr get_raw_nifti_sptr() const; - /// Get image as nifti - std::shared_ptr get_raw_nifti_sptr(); + /// Get image as nifti + std::shared_ptr get_raw_nifti_sptr(); - /// Save to file. Templated so the user can choose the datatype they save to. This defaults - /// to -1, which is the original datatype of that image (stored as _original_datatype). - virtual void write(const std::string &filename, const int datatype) const; + /// Save to file. Templated so the user can choose the datatype they save to. This defaults + /// to -1, which is the original datatype of that image (stored as _original_datatype). + virtual void write(const std::string& filename, const int datatype) const; - /// Write - virtual void write(const std::string &filename) const { this->write(filename,-1); } + /// Write + virtual void write(const std::string& filename) const { this->write(filename, -1); } - /// Get max - float get_max() const; + /// Get max + float get_max() const; - /// Get min - float get_min() const; + /// Get min + float get_min() const; - /// Get mean - float get_mean() const; + /// Get mean + float get_mean() const; - /// Get variance - float get_variance() const; + /// Get variance + float get_variance() const; - /// Get standard deviation - float get_standard_deviation() const; + /// Get standard deviation + float get_standard_deviation() const; - /// Get sum - float get_sum() const; + /// Get sum + float get_sum() const; - /// Get nan count - unsigned get_nan_count() const; + /// Get nan count + unsigned get_nan_count() const; - /// Fill - void fill(const float v); + /// Fill + void fill(const float v); - /// Fill from array - void fill(const dataType *v); + /// Fill from array + void fill(const dataType* v); - /// Fill from array - void fill(const NiftiImageData &im); + /// Fill from array + void fill(const NiftiImageData& im); - /// Get norm - float get_norm(const NiftiImageData&) const; + /// Get norm + float get_norm(const NiftiImageData&) const; - /// Get data size in each dimension - const int* get_dimensions() const; + /// Get data size in each dimension + const int* get_dimensions() const; - /// Get total number of voxels - size_t get_num_voxels() const; + /// Get total number of voxels + size_t get_num_voxels() const; - /// Print header info - void print_header() const; + /// Print header info + void print_header() const; - /// Print multiple header info - static void print_headers(const std::vector &ims); + /// Print multiple header info + static void print_headers(const std::vector& ims); - /// Crop. Set to -1 to leave unchanged - void crop(const int min_index[7], const int max_index[7]); + /// Crop. Set to -1 to leave unchanged + void crop(const int min_index[7], const int max_index[7]); - /// Pad image with value. Give number of voxels to increase in min and max directions. Set values to -1 to leave unchanged - void pad(const int min_index[7], const int max_index[7], const dataType val = 0); + /// Pad image with value. Give number of voxels to increase in min and max directions. Set values to -1 to leave unchanged + void pad(const int min_index[7], const int max_index[7], const dataType val = 0); - /// get 1D index from ND index - int get_1D_index(const int idx[7]) const; + /// get 1D index from ND index + int get_1D_index(const int idx[7]) const; - /// Get original datatype - int get_original_datatype() const { return _original_datatype; } + /// Get original datatype + int get_original_datatype() const { return _original_datatype; } - /// Check if the norms of two images are equal to a given accuracy. - static bool are_equal_to_given_accuracy(const NiftiImageData &im1, const NiftiImageData &im2, const float required_accuracy_compared_to_max); + /// Check if the norms of two images are equal to a given accuracy. + static bool are_equal_to_given_accuracy(const NiftiImageData& im1, const NiftiImageData& im2, const float required_accuracy_compared_to_max); - /// Point is in bounds? - bool is_in_bounds(const int index[7]) const; + /// Point is in bounds? + bool is_in_bounds(const int index[7]) const; - /// Point is in bounds? - bool is_in_bounds(const int index) const; + /// Point is in bounds? + bool is_in_bounds(const int index) const; - /// Images are same size - bool is_same_size(const NiftiImageData &im) const; + /// Images are same size + bool is_same_size(const NiftiImageData& im) const; - /// Do nifti image metadatas match? - static bool do_nifti_image_metadata_match(const NiftiImageData &im1, const NiftiImageData &im2, bool verbose); + /// Do nifti image metadatas match? + static bool do_nifti_image_metadata_match(const NiftiImageData& im1, const NiftiImageData& im2, bool verbose); - /// Dump info of multiple nifti images - static void dump_headers(const std::vector &ims); + /// Dump info of multiple nifti images + static void dump_headers(const std::vector& ims); - /// Dump nifti element - template - static void dump_nifti_element(const std::vector &ims, const std::string &name, const T &call_back); + /// Dump nifti element + template + static void dump_nifti_element(const std::vector& ims, const std::string& name, const T& call_back); - /// Dump nifti element - template - static void dump_nifti_element(const std::vector &ims, const std::string &name, const T &call_back, const unsigned num_elems); + /// Dump nifti element + template + static void dump_nifti_element(const std::vector& ims, const std::string& name, const T& call_back, const unsigned num_elems); - static std::string get_headers(const std::vector*> &ims); - template - static std::string get_nifti_element(const std::vector &ims, const std::string &name, const T &call_back); - template - static std::string get_nifti_element(const std::vector &ims, const std::string &name, const T &call_back, const unsigned num_elems); + static std::string get_headers(const std::vector*>& ims); + template + static std::string get_nifti_element(const std::vector& ims, const std::string& name, const T& call_back); + template + static std::string get_nifti_element(const std::vector& ims, const std::string& name, const T& call_back, const unsigned num_elems); - /// Set the voxel spacing. Requires resampling image, and so interpolation order is required. - /// As per NiftyReg, interpolation_order can be either 0, 1 or 3 meaning nearest neighbor, linear or cubic spline interpolation. - void set_voxel_spacing(const float factors[3], const int interpolation_order); + /// Set the voxel spacing. Requires resampling image, and so interpolation order is required. + /// As per NiftyReg, interpolation_order can be either 0, 1 or 3 meaning nearest neighbor, linear or cubic spline interpolation. + void set_voxel_spacing(const float factors[3], const int interpolation_order); - /// Kernel convolution - void kernel_convolution(const float sigma, NREG_CONV_KERNEL_TYPE conv_type = GAUSSIAN_KERNEL); + /// Kernel convolution + void kernel_convolution(const float sigma, NREG_CONV_KERNEL_TYPE conv_type = GAUSSIAN_KERNEL); - /// Does the image contain any NaNs? - bool get_contains_nans() const { return (this->get_nan_count() > 0); } + /// Does the image contain any NaNs? + bool get_contains_nans() const { return (this->get_nan_count() > 0); } - /// Flip the image along a given axis (Rotation of 180 degrees about axis) - void flip_along_axis(const unsigned axis); + /// Flip the image along a given axis (Rotation of 180 degrees about axis) + void flip_along_axis(const unsigned axis); - /// Mirror the image along a given axis (This will change handedness of image) - void mirror_along_axis(const unsigned axis); + /// Mirror the image along a given axis (This will change handedness of image) + void mirror_along_axis(const unsigned axis); - /// Inner product of two images. - dataType get_inner_product(const NiftiImageData &other) const; + /// Inner product of two images. + dataType get_inner_product(const NiftiImageData& other) const; - /// Normalise image between 0 and 1 - void normalise_zero_and_one(); + /// Normalise image between 0 and 1 + void normalise_zero_and_one(); - /// Standardise (subtract mean and divide by standard deviation). - void standardise(); + /// Standardise (subtract mean and divide by standard deviation). + void standardise(); -protected: + protected: - enum NiftiImageDataType { _general, _3D, _3DTensor, _3DDisp, _3DDef}; + enum NiftiImageDataType { _general, _3D, _3DTensor, _3DDisp, _3DDef }; - enum MathsType { ADD, sub, mul, div}; + enum MathsType { ADD, sub, mul, div }; - /// Image data as a nifti object - std::shared_ptr _nifti_image; + /// Image data as a nifti object + std::shared_ptr _nifti_image; - /// Data - float *_data = NULL; + /// Data + float* _data = NULL; - /// Original datatype - int _original_datatype = -1; + /// Original datatype + int _original_datatype = -1; - /// Check dimensions. Don't require anything for this class. - void check_dimensions(const enum NiftiImageDataType image_type = _general); + /// Check dimensions. Don't require anything for this class. + void check_dimensions(const enum NiftiImageDataType image_type = _general); - /// Set up datatype. Set to float if not already, store the original type. - void set_up_data(const int original_datatype); + /// Set up datatype. Set to float if not already, store the original type. + void set_up_data(const int original_datatype); - /// Add, subract image from another - void maths(const NiftiImageData& c, const MathsType type); + /// Add, subract image from another + void maths(const NiftiImageData& c, const MathsType type); - /// Add, subract, multiply value to image - void maths(const float val, const MathsType type); + /// Add, subract, multiply value to image + void maths(const float val, const MathsType type); - /// Open nifti image - static void open_nifti_image(std::shared_ptr &image, const std::string &filename); + /// Open nifti image + static void open_nifti_image(std::shared_ptr& image, const std::string& filename); - /// Copy nifti image - static void copy_nifti_image(std::shared_ptr &output_image_sptr, const std::shared_ptr &image_to_copy_sptr); + /// Copy nifti image + static void copy_nifti_image(std::shared_ptr& output_image_sptr, const std::shared_ptr& image_to_copy_sptr); -private: + private: - /// Change image datatype with int. Values can be found in nifti1.h (e.g., #define DT_BINARY 1) - void change_datatype(const int datatype); + /// Change image datatype with int. Values can be found in nifti1.h (e.g., #define DT_BINARY 1) + void change_datatype(const int datatype); - /// Change datatype. Templated for desired type. Figures out what current type is then calls doubley templated function below. - template - static void change_datatype(NiftiImageData &im) - { - if (im.get_raw_nifti_sptr()->datatype == DT_BINARY) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_INT8) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_INT16) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_INT32) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT32) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT64) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_UINT8) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_UINT16) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_UINT32) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_INT64) return change_datatype (im); - if (im.get_raw_nifti_sptr()->datatype == DT_UINT64) return change_datatype(im); - if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT128) return change_datatype (im); - - std::stringstream ss; - ss << "change_datatype not implemented for your data type: "; - ss << nifti_datatype_string(im.get_raw_nifti_sptr()->datatype); - ss << " (bytes per voxel: "; - ss << im.get_raw_nifti_sptr()->nbyper << ")."; - throw std::runtime_error(ss.str()); - } - - /// Convert type (performs deep copy) - template - static void change_datatype(NiftiImageData &image) - { - // If the two types are equal, nothing to be done. - if (typeid (newType) == typeid(oldType)) - return; + /// Change datatype. Templated for desired type. Figures out what current type is then calls doubley templated function below. + template + static void change_datatype(NiftiImageData& im) + { + if (im.get_raw_nifti_sptr()->datatype == DT_BINARY) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_INT8) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_INT16) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_INT32) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT32) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT64) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_UINT8) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_UINT16) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_UINT32) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_INT64) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_UINT64) return change_datatype(im); + if (im.get_raw_nifti_sptr()->datatype == DT_FLOAT128) return change_datatype(im); - nifti_image *im = image.get_raw_nifti_sptr().get(); - - // Copy the original array - oldType *originalArray = static_cast(malloc(im->nvox*im->nbyper)); - memcpy(originalArray, im->data, im->nvox*im->nbyper); - // Reset image array - free(im->data); - - // Set the datatype - if (typeid(newType) == typeid(bool)) im->datatype = DT_BINARY; - else if (typeid(newType) == typeid(signed char)) im->datatype = DT_INT8; - else if (typeid(newType) == typeid(signed short)) im->datatype = DT_INT16; - else if (typeid(newType) == typeid(signed int)) im->datatype = DT_INT32; - else if (typeid(newType) == typeid(float)) im->datatype = DT_FLOAT32; - else if (typeid(newType) == typeid(double)) im->datatype = DT_FLOAT64; - else if (typeid(newType) == typeid(unsigned char)) im->datatype = DT_UINT8; - else if (typeid(newType) == typeid(unsigned short)) im->datatype = DT_UINT16; - else if (typeid(newType) == typeid(unsigned int)) im->datatype = DT_UINT32; - else if (typeid(newType) == typeid(signed long long)) im->datatype = DT_INT64; - else if (typeid(newType) == typeid(unsigned long long)) im->datatype = DT_UINT64; - else if (typeid(newType) == typeid(long double)) im->datatype = DT_FLOAT128; - else { std::stringstream ss; - ss << "change_datatype not implemented for your new data type: "; - ss << typeid (newType).name(); + ss << "change_datatype not implemented for your data type: "; + ss << nifti_datatype_string(im.get_raw_nifti_sptr()->datatype); ss << " (bytes per voxel: "; - ss << sizeof(newType) << ")."; + ss << im.get_raw_nifti_sptr()->nbyper << ")."; throw std::runtime_error(ss.str()); } - // Set the nbyper and swap size from the datatype - nifti_datatype_sizes(im->datatype, &im->nbyper, &im->swapsize); - - // Copy data - im->data = static_cast(calloc(im->nvox,sizeof(newType))); - newType *dataPtr = static_cast(im->data); - for (size_t i = 0; i < im->nvox; i++) - dataPtr[i] = newType(originalArray[i]); - - free(originalArray); - return; - } - - // ------------------------------------------------------------------------------ // - // Pure virtual methods from ImageData - // ------------------------------------------------------------------------------ // -public: - /// Clone and return as unique pointer. - std::unique_ptr clone() const - { - return std::unique_ptr(this->clone_impl()); - } - virtual Iterator& begin() - { - _begin.reset(new Iterator(_data)); - return *_begin; - } - virtual Iterator_const& begin() const - { - _begin_const.reset(new Iterator_const(_data)); - return *_begin_const; - } - virtual Iterator& end() - { - _end.reset(new Iterator(_data+_nifti_image->nvox)); - return *_end; - } - virtual Iterator_const& end() const - { - _end_const.reset(new Iterator_const(_data+_nifti_image->nvox)); - return *_end_const; - } -protected: - /// Clone helper function. Don't use. - virtual NiftiImageData* clone_impl() const - { - return new NiftiImageData(*this); - } - virtual ObjectHandle* new_data_container_handle() const - { - return new ObjectHandle - (std::shared_ptr(new NiftiImageData)); - } -public: - unsigned int items() const { return 1; } - /// below all void* are actually float* - virtual void sum (void* ptr) const; - virtual void max (void* ptr) const; - virtual void dot (const DataContainer& a_x, void* ptr) const; - virtual void axpby (const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y); - virtual void xapyb (const DataContainer& a_x, const void* ptr_a, const DataContainer& a_y, const void* ptr_b); - virtual void xapyb (const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b); - virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, - const DataContainer& a_y, const DataContainer& a_b); - virtual float norm() const; - virtual void multiply (const DataContainer& a_x, const DataContainer& a_y); - virtual void divide (const DataContainer& a_x, const DataContainer& a_y); - virtual void maximum(const DataContainer& x, const DataContainer& y); - virtual void minimum(const DataContainer& x, const DataContainer& y); - virtual void power(const DataContainer& x, const DataContainer& y); - virtual void multiply(const DataContainer& a_x, const void* a_y); - virtual void add(const DataContainer& a_x, const void* a_y); - virtual void maximum(const DataContainer& x, const void* a_y); - virtual void minimum(const DataContainer& x, const void* a_y); - virtual void power(const DataContainer& x, const void* a_y); - virtual void exp(const DataContainer& x); - virtual void log(const DataContainer& x); - virtual void sqrt(const DataContainer& x); - virtual void sign(const DataContainer& x); - virtual void abs(const DataContainer& x); - - virtual Dimensions dimensions() const - { - Dimensions dim; - int *d = _nifti_image->dim; - dim["x"] = d[1]; - dim["y"] = d[2]; - dim["z"] = d[3]; - dim["t"] = d[4]; - dim["u"] = d[5]; - dim["v"] = d[6]; - dim["w"] = d[7]; - return dim; - } - void unary_op(const DataContainer& a_x, dataType(*f)(dataType)); - void semibinary_op(const DataContainer& a_x, const void* a_y, dataType(*f)(dataType, dataType)); - void binary_op(const DataContainer& a_x, const DataContainer& a_y, dataType(*f)(dataType, dataType)); - /// Set up the geometrical info. Use qform preferentially over sform. - virtual void set_up_geom_info(); -protected: - mutable std::shared_ptr _begin; - mutable std::shared_ptr _end; - mutable std::shared_ptr _begin_const; - mutable std::shared_ptr _end_const; -}; + /// Convert type (performs deep copy) + template + static void change_datatype(NiftiImageData& image) + { + // If the two types are equal, nothing to be done. + if (typeid (newType) == typeid(oldType)) + return; + + nifti_image* im = image.get_raw_nifti_sptr().get(); + + // Copy the original array + oldType* originalArray = static_cast(malloc(im->nvox * im->nbyper)); + memcpy(originalArray, im->data, im->nvox * im->nbyper); + // Reset image array + free(im->data); + + // Set the datatype + if (typeid(newType) == typeid(bool)) im->datatype = DT_BINARY; + else if (typeid(newType) == typeid(signed char)) im->datatype = DT_INT8; + else if (typeid(newType) == typeid(signed short)) im->datatype = DT_INT16; + else if (typeid(newType) == typeid(signed int)) im->datatype = DT_INT32; + else if (typeid(newType) == typeid(float)) im->datatype = DT_FLOAT32; + else if (typeid(newType) == typeid(double)) im->datatype = DT_FLOAT64; + else if (typeid(newType) == typeid(unsigned char)) im->datatype = DT_UINT8; + else if (typeid(newType) == typeid(unsigned short)) im->datatype = DT_UINT16; + else if (typeid(newType) == typeid(unsigned int)) im->datatype = DT_UINT32; + else if (typeid(newType) == typeid(signed long long)) im->datatype = DT_INT64; + else if (typeid(newType) == typeid(unsigned long long)) im->datatype = DT_UINT64; + else if (typeid(newType) == typeid(long double)) im->datatype = DT_FLOAT128; + else { + std::stringstream ss; + ss << "change_datatype not implemented for your new data type: "; + ss << typeid (newType).name(); + ss << " (bytes per voxel: "; + ss << sizeof(newType) << ")."; + throw std::runtime_error(ss.str()); + } + + // Set the nbyper and swap size from the datatype + nifti_datatype_sizes(im->datatype, &im->nbyper, &im->swapsize); + + // Copy data + im->data = static_cast(calloc(im->nvox, sizeof(newType))); + newType* dataPtr = static_cast(im->data); + for (size_t i = 0; i < im->nvox; i++) + dataPtr[i] = newType(originalArray[i]); + + free(originalArray); + return; + } + + // ------------------------------------------------------------------------------ // + // Pure virtual methods from ImageData + // ------------------------------------------------------------------------------ // + public: + /// Clone and return as unique pointer. + std::unique_ptr clone() const + { + return std::unique_ptr(this->clone_impl()); + } + virtual Iterator& begin() + { + _begin.reset(new Iterator(_data)); + return *_begin; + } + virtual Iterator_const& begin() const + { + _begin_const.reset(new Iterator_const(_data)); + return *_begin_const; + } + virtual Iterator& end() + { + _end.reset(new Iterator(_data + _nifti_image->nvox)); + return *_end; + } + virtual Iterator_const& end() const + { + _end_const.reset(new Iterator_const(_data + _nifti_image->nvox)); + return *_end_const; + } + protected: + /// Clone helper function. Don't use. + virtual NiftiImageData* clone_impl() const + { + return new NiftiImageData(*this); + } + virtual ObjectHandle* new_data_container_handle() const + { + return new ObjectHandle + (std::shared_ptr(new NiftiImageData)); + } + public: + unsigned int items() const { return 1; } + virtual float norm() const; + virtual float sum() const; + virtual float max() const; + virtual float dot(const DataContainer& a_x) const; + virtual void axpby(dataType a, const DataContainer& a_x, dataType b, const DataContainer& a_y); + virtual void xapyb(const DataContainer& a_x, dataType a, const DataContainer& a_y, dataType b); + virtual void xapyb(const DataContainer& a_x, dataType a, const DataContainer& a_y, const DataContainer& a_b); + virtual void xapyb(const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b); + + virtual Dimensions dimensions() const + { + Dimensions dim; + int* d = _nifti_image->dim; + dim["x"] = d[1]; + dim["y"] = d[2]; + dim["z"] = d[3]; + dim["t"] = d[4]; + dim["u"] = d[5]; + dim["v"] = d[6]; + dim["w"] = d[7]; + return dim; + } + virtual void unary_op(const DataContainer& a_x, dataType(*f)(dataType)); + virtual void semibinary_op(const DataContainer& a_x, dataType y, dataType(*f)(dataType, dataType)); + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, dataType(*f)(dataType, dataType)); + /// Set up the geometrical info. Use qform preferentially over sform. + virtual void set_up_geom_info(); + protected: + mutable std::shared_ptr _begin; + mutable std::shared_ptr _end; + mutable std::shared_ptr _begin_const; + mutable std::shared_ptr _end_const; + }; } diff --git a/src/Registration/cReg/include/sirf/Reg/NiftiImageData3D.h b/src/Registration/cReg/include/sirf/Reg/NiftiImageData3D.h index 55ef6fb1c..3e45df09e 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftiImageData3D.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftiImageData3D.h @@ -68,7 +68,7 @@ class NiftiImageData3D : public NiftiImageData : NiftiImageData(data, geom) { this->check_dimensions(this->_3D); } /// Construct from any other image data (e.g., STIRImageData) - NiftiImageData3D(const ImageData& id) + NiftiImageData3D(const ImageData& id) : NiftiImageData(id) { this->check_dimensions(this->_3D); } virtual ObjectHandle* new_data_container_handle() const diff --git a/src/Registration/cReg/include/sirf/Reg/NiftiImageData3DTensor.h b/src/Registration/cReg/include/sirf/Reg/NiftiImageData3DTensor.h index d7862cb21..bdf48b1a0 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftiImageData3DTensor.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftiImageData3DTensor.h @@ -91,13 +91,13 @@ class NiftiImageData3DTensor : public NiftiImageData void flip_component(const int dim); /// Tensor component maths - void tensor_component_maths(const int dim, const std::shared_ptr &scalar_im_sptr, const typename NiftiImageData::MathsType maths_type); + void tensor_component_maths(const int dim, const std::shared_ptr > &scalar_im_sptr, const typename NiftiImageData::MathsType maths_type); /// Multiply tensor component by image - void multiply_tensor_component(const int dim, const std::shared_ptr &scalar_im_sptr); + void multiply_tensor_component(const int dim, const std::shared_ptr > &scalar_im_sptr); /// Add image to tensor component - void add_to_tensor_component(const int dim, const std::shared_ptr &scalar_im_sptr); + void add_to_tensor_component(const int dim, const std::shared_ptr > &scalar_im_sptr); virtual ObjectHandle* new_data_container_handle() const { diff --git a/src/Registration/cReg/include/sirf/Reg/NiftyRegistration.h b/src/Registration/cReg/include/sirf/Reg/NiftyRegistration.h index 129852693..7ccae8db8 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftyRegistration.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftyRegistration.h @@ -64,10 +64,10 @@ class NiftyRegistration : public NiftiBasedRegistration void set_parameter(const std::string &par, const std::string &arg1 = "", const std::string &arg2 = ""); /// Set reference mask - void set_reference_mask(const std::shared_ptr reference_mask_sptr) { _reference_mask_sptr = reference_mask_sptr; } + void set_reference_mask(const std::shared_ptr > reference_mask_sptr) { _reference_mask_sptr = reference_mask_sptr; } /// Set floating mask - void set_floating_mask(const std::shared_ptr floating_mask_sptr) { _floating_mask_sptr = floating_mask_sptr; } + void set_floating_mask(const std::shared_ptr > floating_mask_sptr) { _floating_mask_sptr = floating_mask_sptr; } protected: @@ -87,9 +87,9 @@ class NiftyRegistration : public NiftiBasedRegistration std::string _parameter_filename; /// Floating mask - std::shared_ptr _floating_mask_sptr; + std::shared_ptr > _floating_mask_sptr; /// Reference mask - std::shared_ptr _reference_mask_sptr; + std::shared_ptr > _reference_mask_sptr; /// Floating mask (as NiftiImageData3D) std::shared_ptr > _floating_mask_nifti_sptr; diff --git a/src/Registration/cReg/include/sirf/Reg/NiftyResampler.h b/src/Registration/cReg/include/sirf/Reg/NiftyResampler.h index 5136330f0..50b9eaeb4 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftyResampler.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftyResampler.h @@ -112,16 +112,16 @@ class NiftyResampler : public Resampler virtual void process(); /// Do the forward transformation - virtual std::shared_ptr forward(const std::shared_ptr input_sptr); + virtual std::shared_ptr forward(const std::shared_ptr input_sptr); /// Do the forward transformation - virtual void forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); + virtual void forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); /// Do the adjoint transformation - virtual std::shared_ptr adjoint(const std::shared_ptr input_sptr); + virtual std::shared_ptr adjoint(const std::shared_ptr input_sptr); /// Do the adjoint transformation - virtual void adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); + virtual void adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); protected: diff --git a/src/Registration/cReg/include/sirf/Reg/Registration.h b/src/Registration/cReg/include/sirf/Reg/Registration.h index 99ee62e02..268a6eb01 100644 --- a/src/Registration/cReg/include/sirf/Reg/Registration.h +++ b/src/Registration/cReg/include/sirf/Reg/Registration.h @@ -37,7 +37,7 @@ namespace sirf { /// Forward declarations template class Transformation; -class ImageData; +template class ImageData; /*! \ingroup Registration @@ -79,13 +79,13 @@ class Registration virtual ~Registration() {} /// Set reference image - void set_reference_image(const std::shared_ptr reference_image_sptr); + void set_reference_image(const std::shared_ptr > reference_image_sptr); /// Set floating image. Will clear any previous floating images. - void set_floating_image(const std::shared_ptr floating_image_sptr); + void set_floating_image(const std::shared_ptr > floating_image_sptr); /// Add floating image - void add_floating_image(const std::shared_ptr floating_image_sptr); + void add_floating_image(const std::shared_ptr > floating_image_sptr); /// Set reference image filename. Will be read as NiftiImageData. void set_reference_image_filename(const std::string &filename); @@ -103,7 +103,7 @@ class Registration virtual void process() = 0; /// Get registered image - virtual const std::shared_ptr get_output_sptr(const unsigned idx = 0) const { return _warped_images.at(idx); } + virtual const std::shared_ptr > get_output_sptr(const unsigned idx = 0) const { return _warped_images.at(idx); } /// Get forward deformation field image virtual const std::shared_ptr > get_deformation_field_forward_sptr(const unsigned idx = 0) const { return _def_fwd_images.at(idx); } @@ -123,11 +123,11 @@ class Registration virtual void check_parameters() const; /// Reference image - std::shared_ptr _reference_image_sptr; + std::shared_ptr > _reference_image_sptr; /// Floating image - std::vector > _floating_images; + std::vector > > _floating_images; /// Warped image - std::vector > _warped_images; + std::vector > > _warped_images; /// Forward deformation field image std::vector > > _def_fwd_images; diff --git a/src/Registration/cReg/include/sirf/Reg/Resample.h b/src/Registration/cReg/include/sirf/Reg/Resample.h index 032400478..084df09a4 100644 --- a/src/Registration/cReg/include/sirf/Reg/Resample.h +++ b/src/Registration/cReg/include/sirf/Reg/Resample.h @@ -38,7 +38,8 @@ namespace sirf { // Forward declarations template class Transformation; -class ImageData; +//template class ImageData; +class DataContainer; /*! \file @@ -72,10 +73,10 @@ class Resampler virtual ~Resampler() {} /// Set reference image. This is the image that would be the reference if you were doing a forward transformation. - virtual void set_reference_image(const std::shared_ptr reference_image_sptr); + virtual void set_reference_image(const std::shared_ptr reference_image_sptr); /// Set floating image. This is the image that would be the floating if you were doing a forward transformation. - virtual void set_floating_image(const std::shared_ptr floating_image_sptr); + virtual void set_floating_image(const std::shared_ptr floating_image_sptr); /// Add transformation virtual void add_transformation(const std::shared_ptr > transformation_sptr); @@ -108,25 +109,25 @@ class Resampler virtual void process() = 0; /// Get output - const std::shared_ptr get_output_sptr() const { return _output_image_sptr; } + const std::shared_ptr get_output_sptr() const { return _output_image_sptr; } /// Do the forward transformation - virtual std::shared_ptr forward(const std::shared_ptr input_sptr) = 0; + virtual std::shared_ptr forward(const std::shared_ptr input_sptr) = 0; /// Do the forward transformation - virtual void forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) = 0; + virtual void forward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) = 0; /// Do the adjoint transformation - virtual std::shared_ptr adjoint(const std::shared_ptr input_sptr) = 0; + virtual std::shared_ptr adjoint(const std::shared_ptr input_sptr) = 0; /// Do the adjoint transformation - virtual void adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) = 0; + virtual void adjoint(std::shared_ptr output_sptr, const std::shared_ptr input_sptr) = 0; /// Backward. Alias for Adjoint - virtual std::shared_ptr backward(const std::shared_ptr input_sptr); + virtual std::shared_ptr backward(const std::shared_ptr input_sptr); /// Backward. Alias for Adjoint - virtual void backward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); + virtual void backward(std::shared_ptr output_sptr, const std::shared_ptr input_sptr); protected: @@ -143,9 +144,9 @@ class Resampler virtual void check_parameters(); /// Reference image - std::shared_ptr _reference_image_sptr; + std::shared_ptr _reference_image_sptr; /// Floating image - std::shared_ptr _floating_image_sptr; + std::shared_ptr _floating_image_sptr; /// Transformations (could be mixture of affine, displacements, deformations). std::vector > > _transformations; @@ -154,7 +155,7 @@ class Resampler InterpolationType _interpolation_type; /// Output image - std::shared_ptr _output_image_sptr; + std::shared_ptr _output_image_sptr; /// Padding value float _padding_value = 0; diff --git a/src/Synergistic/cSyn/include/sirf/Syn/utilities.h b/src/Synergistic/cSyn/include/sirf/Syn/utilities.h index 9e1b0a781..e727ed7ce 100644 --- a/src/Synergistic/cSyn/include/sirf/Syn/utilities.h +++ b/src/Synergistic/cSyn/include/sirf/Syn/utilities.h @@ -34,23 +34,23 @@ namespace sirf { class ImageDataWrap { public: ImageDataWrap(const std::string &filename, const std::string &engine, bool verbose); - ImageData& data() + DataContainer& data() { return *img_sptr_; } - const ImageData& data() const + const DataContainer& data() const { return *img_sptr_; } - std::shared_ptr data_sptr() + std::shared_ptr data_sptr() { return img_sptr_; } - const std::shared_ptr data_sptr() const + const std::shared_ptr data_sptr() const { return img_sptr_; } private: - std::shared_ptr img_sptr_; + std::shared_ptr img_sptr_; }; } diff --git a/src/Synergistic/cSyn/tests/test_conv_img.cpp b/src/Synergistic/cSyn/tests/test_conv_img.cpp index 5183f1f6c..b04f5daaa 100644 --- a/src/Synergistic/cSyn/tests/test_conv_img.cpp +++ b/src/Synergistic/cSyn/tests/test_conv_img.cpp @@ -42,7 +42,7 @@ int main(int argc, char* argv[]) ImageDataWrap imw(filename, eng_in, true); std::cout << "converting " << eng_in.c_str() << " image to " << eng_out.c_str() << " image...\n"; - const ImageData& im_in = imw.data(); + const ImageData& im_in = (const ImageData&)imw.data(); if (eng_out == std::string("Reg")) { NiftiImageData3D im_out(im_in); if (im_out == im_in) { diff --git a/src/Synergistic/cSyn/utilities.cpp b/src/Synergistic/cSyn/utilities.cpp index e982b8163..e05ac0bb7 100644 --- a/src/Synergistic/cSyn/utilities.cpp +++ b/src/Synergistic/cSyn/utilities.cpp @@ -38,28 +38,34 @@ using namespace sirf; ImageDataWrap::ImageDataWrap(const std::string &filename, const std::string &engine, bool verbose) { + std::shared_ptr gi_sptr; if (strcmp(engine.c_str(), "Reg") == 0) { std::shared_ptr > nifti_sptr = std::make_shared >(filename); if (verbose) nifti_sptr->print_header(); img_sptr_ = nifti_sptr; + gi_sptr = nifti_sptr->get_geom_info_sptr(); } else if (strcmp(engine.c_str(), "STIR") == 0) { - img_sptr_ = std::make_shared(filename); + std::shared_ptr stir_sptr = std::make_shared(filename); + img_sptr_ = stir_sptr; + gi_sptr = stir_sptr->get_geom_info_sptr(); + //img_sptr_ = std::make_shared(filename); } else if (strcmp(engine.c_str(), "Gadgetron") == 0) { std::shared_ptr gadgetron_sptr(new GadgetronImagesVector); gadgetron_sptr->read(filename); if (verbose) gadgetron_sptr->print_header(0); img_sptr_ = gadgetron_sptr; + gi_sptr = gadgetron_sptr->get_geom_info_sptr(); } else throw std::runtime_error("unknown engine - " + engine + ".\n"); // If verbose print geom info if (verbose) { - std::shared_ptr gi_sptr = - img_sptr_->get_geom_info_sptr(); + //std::shared_ptr gi_sptr = + // img_sptr_->get_geom_info_sptr(); if (gi_sptr.get()) gi_sptr->print_info(); } diff --git a/src/Synergistic/sirf_convert_image_type.cpp b/src/Synergistic/sirf_convert_image_type.cpp index 8e639e163..d6b25d272 100644 --- a/src/Synergistic/sirf_convert_image_type.cpp +++ b/src/Synergistic/sirf_convert_image_type.cpp @@ -33,13 +33,14 @@ limitations under the License. using namespace sirf; -static const std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine, const bool verbose) +static const std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine, const bool verbose) { - std::shared_ptr img_sptr; + std::shared_ptr img_sptr; if (strcmp(engine.c_str(), "Reg") == 0) { std::shared_ptr > nifti_sptr = std::make_shared >(filename); if (verbose) nifti_sptr->print_header(); + if (verbose) nifti_sptr->get_geom_info_sptr()->print_info(); img_sptr = nifti_sptr; } else if (strcmp(engine.c_str(), "STIR") == 0) { @@ -50,22 +51,25 @@ static const std::shared_ptr image_as_sptr(const std::string &filenam std::shared_ptr gadgetron_sptr(new GadgetronImagesVector); gadgetron_sptr->read(filename); if (verbose) gadgetron_sptr->print_header(0); + if (verbose) gadgetron_sptr->get_geom_info_sptr()->print_info(); img_sptr = gadgetron_sptr; } else throw std::runtime_error("unknown engine - " + engine + ".\n"); - // If verbose print geom info - if (verbose) img_sptr->get_geom_info_sptr()->print_info(); + //// If verbose print geom info + //if (verbose) img_sptr->get_geom_info_sptr()->print_info(); // return return img_sptr; } -static void convert_and_write_image(const std::string &filename, const std::string &engine, const std::shared_ptr &in_img_sptr, const std::string ¶m_file, const bool verbose) +static void convert_and_write_image(const std::string &filename, const std::string &engine, const std::shared_ptr &in_img_sptr, const std::string ¶m_file, const bool verbose) { if (strcmp(engine.c_str(), "Reg") == 0) { - NiftiImageData im(*in_img_sptr); + auto im_sptr = std::dynamic_pointer_cast >(in_img_sptr); + auto im = *im_sptr; + //NiftiImageData im(*in_img_sptr); im.write(filename); if (verbose) { im.print_header(); @@ -73,7 +77,9 @@ static void convert_and_write_image(const std::string &filename, const std::stri } } else if (strcmp(engine.c_str(), "STIR") == 0) { - STIRImageData im(*in_img_sptr); + const STIRImageData* ptr = dynamic_cast(in_img_sptr.get()); + STIRImageData im(*ptr); + // STIRImageData im(*in_img_sptr); if (param_file.empty()) im.write(filename); else @@ -174,7 +180,7 @@ int main(int argc, char* argv[]) } // Read input image - std::shared_ptr in_img_sptr = image_as_sptr(in_filename, in_engine, verbose); + std::shared_ptr in_img_sptr = image_as_sptr(in_filename, in_engine, verbose); // Convert and write convert_and_write_image(out_filename, out_engine, in_img_sptr, param_file, verbose); diff --git a/src/Synergistic/sirf_registration.cpp b/src/Synergistic/sirf_registration.cpp index 1d38312b4..10be94d84 100644 --- a/src/Synergistic/sirf_registration.cpp +++ b/src/Synergistic/sirf_registration.cpp @@ -46,7 +46,7 @@ enum Algorithm { SPM }; -static std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine) +static std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine) { if (strcmp(engine.c_str(), "Reg") == 0) return std::make_shared >(filename); @@ -287,21 +287,21 @@ int main(int argc, char* argv[]) // Ref if (ref_str.empty()) throw std::runtime_error("--ref not set"); - reg->set_reference_image(image_as_sptr(ref_str,ref_eng_str)); + reg->set_reference_image(std::dynamic_pointer_cast >(image_as_sptr(ref_str,ref_eng_str))); // Flo if (flo_strs.empty()) throw std::runtime_error("--flo not set"); for (unsigned i=0; iadd_floating_image(image_as_sptr(flo_strs.at(i).first,flo_strs.at(i).second)); + reg->add_floating_image(std::dynamic_pointer_cast >(image_as_sptr(flo_strs.at(i).first,flo_strs.at(i).second))); // rmask if (!rmask_str.empty()) { if (algo == SPM) throw std::runtime_error("--rmask not available for spm"); - std::dynamic_pointer_cast >(reg)->set_reference_mask(image_as_sptr(rmask_str,rmask_eng_str)); + std::dynamic_pointer_cast >(reg)->set_reference_mask(std::dynamic_pointer_cast >(image_as_sptr(rmask_str,rmask_eng_str))); } // fmask if (!fmask_str.empty()) { if (algo == SPM) throw std::runtime_error("--fmask not available for spm"); - std::dynamic_pointer_cast >(reg)->set_reference_mask(image_as_sptr(fmask_str,fmask_eng_str)); + std::dynamic_pointer_cast >(reg)->set_reference_mask(std::dynamic_pointer_cast >(image_as_sptr(fmask_str,fmask_eng_str))); } // print if (print) { diff --git a/src/Synergistic/sirf_resample.cpp b/src/Synergistic/sirf_resample.cpp index ad79db96c..b86dd135c 100644 --- a/src/Synergistic/sirf_resample.cpp +++ b/src/Synergistic/sirf_resample.cpp @@ -41,7 +41,7 @@ If multiple transformations are given, they will be applied in the order they we using namespace sirf; -static std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine) +static std::shared_ptr image_as_sptr(const std::string &filename, const std::string &engine) { if (strcmp(engine.c_str(), "Nifti") == 0) return std::make_shared >(filename); @@ -217,24 +217,24 @@ int main(int argc, char* argv[]) err("Error: -flo required."); // Get images as NiftiImages - std::shared_ptr ref = image_as_sptr(ref_filename,eng_ref); - std::shared_ptr flo = image_as_sptr(flo_filename,eng_flo); + std::shared_ptr ref = image_as_sptr(ref_filename,eng_ref); + std::shared_ptr flo = image_as_sptr(flo_filename,eng_flo); // Resample std::shared_ptr > res = algo_as_sptr(algo); - res->set_reference_image(ref); - res->set_floating_image(flo); + res->set_reference_image(std::dynamic_pointer_cast>(ref)); + res->set_floating_image(std::dynamic_pointer_cast>(flo)); for (size_t i=0; iadd_transformation(trans[i]); res->set_interpolation_type(interp); if (pad_set) res->set_padding_value(pad); - std::shared_ptr output_sptr; + std::shared_ptr output_sptr; if (forward) - output_sptr = res->forward(flo); + output_sptr = res->forward(std::dynamic_pointer_cast>(flo)); else - output_sptr = res->adjoint(ref); + output_sptr = res->adjoint(std::dynamic_pointer_cast>(ref)); output_sptr->write(output); } diff --git a/src/Synergistic/tests/test_cSynergistic.cpp b/src/Synergistic/tests/test_cSynergistic.cpp index e81f66857..97b133d07 100644 --- a/src/Synergistic/tests/test_cSynergistic.cpp +++ b/src/Synergistic/tests/test_cSynergistic.cpp @@ -154,6 +154,7 @@ int main(int argc, char* argv[]) recon.add_gadget("gadget_" + std::to_string(i), gadgets[i]); recon.process(raw_mr); std::shared_ptr ismrmrd_im_sptr = recon.get_output(); + //std::cout << "ok\n" << std::flush; if (!ismrmrd_im_sptr->is_complex()) throw std::runtime_error("Expected output of reconstruction to be complex"); @@ -166,27 +167,37 @@ int main(int argc, char* argv[]) cmplx_flt.imag(cmplx_flt.real() / 2.f); *iter = NumRef((void *)&cmplx_flt, NumberType::CXFLOAT); } + //std::cout << "ok\n" << std::flush; // Convert the complex image to two niftis std::shared_ptr > real_sptr, imag_sptr; NiftiImageData::construct_NiftiImageData_from_complex_im(real_sptr,imag_sptr,ismrmrd_im_sptr); + //std::cout << imag_sptr->data_type() << '\n'; + //std::cout << "ok1\n" << std::flush; real_sptr->write("results/real"); + //std::cout << "ok2\n"; imag_sptr->write("results/imag"); + //std::cout << "ok3\n"; // Create affine transformation std::shared_ptr > tm_sptr = std::make_shared >(); (*tm_sptr)[0][3] = 2.f; + //std::cout << "ok2\n" << std::flush; // Resample the complex data NiftyResampler res_complex; res_complex.set_reference_image(ismrmrd_im_sptr); res_complex.set_floating_image(ismrmrd_im_sptr); + //std::cout << "ok3\n" << std::flush; res_complex.set_interpolation_type_to_linear(); res_complex.add_transformation(tm_sptr); - std::shared_ptr forward_cplx_sptr = res_complex.forward(ismrmrd_im_sptr); - std::shared_ptr adjoint_cplx_sptr = res_complex.adjoint(ismrmrd_im_sptr); + //std::cout << "ok4\n" << std::flush; + std::shared_ptr forward_cplx_sptr = res_complex.forward(ismrmrd_im_sptr); + //std::cout << "ok5\n" << std::flush; + std::shared_ptr adjoint_cplx_sptr = res_complex.adjoint(ismrmrd_im_sptr); + //std::cout << "ok6\n" << std::flush; // Get the output std::shared_ptr > forward_cplx_real_sptr, forward_cplx_imag_sptr, adjoint_cplx_real_sptr, adjoint_cplx_imag_sptr; @@ -208,6 +219,7 @@ int main(int argc, char* argv[]) std::dynamic_pointer_cast >(res_real.forward(real_sptr)); std::shared_ptr > adjoint_real_sptr = std::dynamic_pointer_cast >(res_real.adjoint(real_sptr)); + //std::cout << "ok7\n" << std::flush; NiftyResampler res_imag; res_imag.set_reference_image(imag_sptr); @@ -218,19 +230,27 @@ int main(int argc, char* argv[]) std::dynamic_pointer_cast >(res_imag.forward(imag_sptr)); std::shared_ptr > adjoint_imag_sptr = std::dynamic_pointer_cast >(res_imag.adjoint(imag_sptr)); + //std::cout << "ok8\n" << std::flush; // Compare that the real and imaginary parts match regardless // of whether they were resampled separately or together. + //std::cout << (*forward_real_sptr == *forward_cplx_real_sptr) << '\n' << std::flush; + //std::cout << (*forward_imag_sptr == *forward_cplx_imag_sptr) << '\n' << std::flush; + //std::cout << (*adjoint_real_sptr == *adjoint_cplx_real_sptr) << '\n' << std::flush; + //std::cout << (*adjoint_imag_sptr == *adjoint_cplx_imag_sptr) << '\n' << std::flush; if (*forward_real_sptr != *forward_cplx_real_sptr || *forward_imag_sptr != *forward_cplx_imag_sptr) throw std::runtime_error("NiftyResampler forward failed for complex data"); + //std::cout << "ok9\n" << std::flush; if (*adjoint_real_sptr != *adjoint_cplx_real_sptr || *adjoint_imag_sptr != *adjoint_cplx_imag_sptr) throw std::runtime_error("NiftyResampler adjoint failed for complex data"); + //std::cout << "ok10\n" << std::flush; std::cout << "// ----------------------------------------------------------------------- //\n"; std::cout << "// Finished complex resampler test.\n"; std::cout << "//------------------------------------------------------------------------ //\n"; + std::cout << std::flush; } // Test MR reorient @@ -242,11 +262,13 @@ int main(int argc, char* argv[]) // Read ISMRMRD image std::shared_ptr G1_sptr = std::make_shared(); G1_sptr->read(mr_recon_h5_filename); + //std::cout << G1_sptr->is_complex() << " ok 1\n" << std::flush; // Convert ISMRMRD image to nifti std::shared_ptr > G1_nii_sptr = std::make_shared >(*G1_sptr); std::shared_ptr raw_nii_sptr = G1_nii_sptr->get_raw_nifti_sptr(); + //std::cout << "ok 2\n" << std::flush; // Affine transformation as translation by integer num voxels (so no interpolation) std::array trans = { G1_sptr->get_geom_info_sptr()->get_spacing()[0] * 2.f, @@ -280,15 +302,26 @@ int main(int argc, char* argv[]) reg_checkAndCorrectDimension(raw_nii_sptr.get()); // Re-set up geom info G1_nii_sptr->set_up_geom_info(); + //std::cout << "ok 3\n" << std::flush; // Reorient gadgetron image with modified G1_nii_sptr's geom info std::shared_ptr G2_sptr = G1_sptr->clone(); + /* + std::shared_ptr x = std::dynamic_pointer_cast(G2_sptr); + std::cout << (size_t)x.get() << '\n' << std::flush; + std::shared_ptr y = std::dynamic_pointer_cast(x); + std::cout << (size_t)y.get() << '\n' << std::flush; + std::shared_ptr > z = std::dynamic_pointer_cast >(x); + std::cout << (size_t)z.get() << '\n' << std::flush; + */ G2_sptr->reorient(*G1_nii_sptr->get_geom_info_sptr()); + //std::cout << "ok 4\n" << std::flush; std::cout << "\n original:\n"; G1_sptr->get_geom_info_sptr()->print_info(); std::cout << "\n resampled:\n"; G2_sptr->get_geom_info_sptr()->print_info(); + //std::cout << "ok 5\n" << std::flush; // Now resampled G2 back to G1 using inverse TM, should be the same NiftyResampler res; @@ -297,15 +330,22 @@ int main(int argc, char* argv[]) res.set_padding_value(0.f); res.set_interpolation_type_to_linear(); res.add_transformation(std::make_shared >(trans_sptr->get_inverse())); - std::shared_ptr resampled_G2_sptr = res.forward(G2_sptr); + //std::cout << "ok 6-\n" << std::flush; + std::shared_ptr resampled_G2_sptr = res.forward(G2_sptr); + //std::cout << "ok 6\n" << std::flush; std::cout << "\n reoriented back to original space:\n"; - resampled_G2_sptr->get_geom_info_sptr()->print_info(); - - if (NiftiImageData(*G1_sptr) != NiftiImageData(*resampled_G2_sptr)) + std::shared_ptr sptr_im2 = std::dynamic_pointer_cast(resampled_G2_sptr); + sptr_im2->get_geom_info_sptr()->print_info(); + //resampled_G2_sptr->get_geom_info_sptr()->print_info(); + //std::cout << "ok 7\n" << std::flush; + + auto sptr_G1 = std::dynamic_pointer_cast >(G1_sptr); + auto sptr_G2 = std::dynamic_pointer_cast >(resampled_G2_sptr); + //if (NiftiImageData(*G1_sptr) != NiftiImageData(*resampled_G2_sptr)) + if (*sptr_G1 != *sptr_G2) throw std::runtime_error("GadgetronImagesVector::reorient test failed"); - std::cout << "// ----------------------------------------------------------------------- //\n"; std::cout << "// Finished GadgetronImageData reorient test.\n"; std::cout << "//------------------------------------------------------------------------ //\n"; diff --git a/src/common/ImageData.cpp b/src/common/ImageData.cpp index d49b217b6..6d22e3e88 100644 --- a/src/common/ImageData.cpp +++ b/src/common/ImageData.cpp @@ -22,19 +22,21 @@ limitations under the License. using namespace sirf; -void ImageData::reorient(const VoxelisedGeometricalInfo3D &) -{ - throw std::runtime_error("ImageData::reorient not yet implemented for your image type."); -} +//template +//void ImageData::reorient(const VoxelisedGeometricalInfo3D & geom_info) +//{ +// throw std::runtime_error("ImageData::reorient not yet implemented for your image type."); +//} -bool ImageData::can_reorient(const VoxelisedGeometricalInfo3D &geom_1, const VoxelisedGeometricalInfo3D &geom_2, const bool throw_error) -{ - // If size and spacing match, return true - if (geom_1.get_size() == geom_2.get_size()) - return true; - // Else (and error desired), print error - if (throw_error) - throw std::runtime_error("ImageData::can_reorient: num voxels do not match."); - // Else, return false - return false; -} \ No newline at end of file +//template +//bool ImageData::can_reorient(const VoxelisedGeometricalInfo3D &geom_1, const VoxelisedGeometricalInfo3D &geom_2, const bool throw_error) +//{ +// // If size and spacing match, return true +// if (geom_1.get_size() == geom_2.get_size()) +// return true; +// // Else (and error desired), print error +// if (throw_error) +// throw std::runtime_error("ImageData::can_reorient: num voxels do not match."); +// // Else, return false +// return false; +//} \ No newline at end of file diff --git a/src/common/SIRF.py b/src/common/SIRF.py index b470b8a86..9e06c3903 100644 --- a/src/common/SIRF.py +++ b/src/common/SIRF.py @@ -286,7 +286,7 @@ def add(self, other, out=None): z.handle = pysirf.cSIRF_axpby(one.ctypes.data, self.handle, one.ctypes.data, other.handle) check_status(z.handle) else: - try_calling(pysirf.cSIRF_axpbyAlt(one.ctypes.data, self.handle, one.ctypes.data, other.handle, z.handle)) + try_calling(pysirf.cSIRF_compute_axpby(one.ctypes.data, self.handle, one.ctypes.data, other.handle, z.handle)) return z def subtract(self, other, out=None): @@ -316,7 +316,7 @@ def subtract(self, other, out=None): else: assert_validities(self, out) z = out - try_calling(pysirf.cSIRF_axpbyAlt \ + try_calling(pysirf.cSIRF_compute_axpby \ (pl_one.ctypes.data, self.handle, mn_one.ctypes.data, other.handle, z.handle)) def multiply(self, other, out=None): @@ -414,13 +414,13 @@ def sapyb(self, a, y, b, out=None, **kwargs): if out is None: z.handle = pysirf.cSIRF_axpby(alpha.ctypes.data, self.handle, beta.ctypes.data, y.handle) else: - try_calling(pysirf.cSIRF_axpbyAlt(alpha.ctypes.data, self.handle, beta.ctypes.data, y.handle, z.handle)) + try_calling(pysirf.cSIRF_compute_axpby(alpha.ctypes.data, self.handle, beta.ctypes.data, y.handle, z.handle)) else: #a is scalar, b is array if out is None: z.handle = pysirf.cSIRF_XapYB(self.handle, alpha.ctypes.data, y.handle, b.handle) else: - try_calling(pysirf.cSIRF_XapYBAlt(self.handle, alpha.ctypes.data, y.handle, b.handle, z.handle)) + try_calling(pysirf.cSIRF_compute_XapYB(self.handle, alpha.ctypes.data, y.handle, b.handle, z.handle)) else: assert_validities(self, a) if isinstance(b, Number): @@ -429,14 +429,14 @@ def sapyb(self, a, y, b, out=None, **kwargs): if out is None: z.handle = pysirf.cSIRF_XapYB(y.handle, beta.ctypes.data, self.handle, a.handle) else: - try_calling(pysirf.cSIRF_XapYBAlt(y.handle, beta.ctypes.data, self.handle, a.handle, z.handle)) + try_calling(pysirf.cSIRF_compute_XapYB(y.handle, beta.ctypes.data, self.handle, a.handle, z.handle)) else: #a is array, b is array assert_validities(self, b) if out is None: z.handle = pysirf.cSIRF_xapyb(self.handle, a.handle, y.handle, b.handle) else: - try_calling(pysirf.cSIRF_xapybAlt(self.handle, a.handle, y.handle, b.handle, z.handle)) + try_calling(pysirf.cSIRF_compute_xapyb(self.handle, a.handle, y.handle, b.handle, z.handle)) if out is None: check_status(z.handle) diff --git a/src/common/csirf.cpp b/src/common/csirf.cpp index 05b80488a..57f918005 100644 --- a/src/common/csirf.cpp +++ b/src/common/csirf.cpp @@ -34,6 +34,17 @@ using namespace sirf; #define NEW_OBJECT_HANDLE(T) new ObjectHandle(shared_ptr(new T)) #define SPTR_FROM_HANDLE(Object, X, H) \ shared_ptr X; getObjectSptrFromHandle(H, X); +#define SELECT_DATA_CASE(DataType, Method, Args, ...) \ + if (DataType == "float") \ + Method(Args, ##__VA_ARGS__); \ + else if (DataType == "complex_float") \ + Method(Args, ##__VA_ARGS__); \ + else { \ + std::string err = "??? Unsupported data type "; \ + err += DataType; \ + throw(err.c_str()); \ + } \ + static void* @@ -141,43 +152,81 @@ cSIRF_norm(const void* ptr_x) CATCH; } +template +void +compute_dot_templ(const void* ptr_x, const void* ptr_y, void* ptr_z) +{ + auto const& x = objectFromHandle >(ptr_x); + auto const& y = objectFromHandle >(ptr_y); + *static_cast(ptr_z) = x.dot(y); +} + extern "C" void* cSIRF_compute_dot(const void* ptr_x, const void* ptr_y, void* ptr_z) { try { - auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); - x.dot(y, ptr_z); + auto const& base = objectFromHandle(ptr_x); + SELECT_DATA_CASE(base.data_type(), compute_dot_templ, ptr_x, ptr_y, ptr_z); return new DataHandle; } CATCH; } +template +void +compute_sum_templ(const void* ptr_x, void* ptr_z) +{ + auto const& x = objectFromHandle >(ptr_x); + *static_cast(ptr_z) = x.sum(); +} + extern "C" void* cSIRF_compute_sum(const void* ptr_x, void* ptr_z) { try { - auto const& x = objectFromHandle(ptr_x); - x.sum(ptr_z); + auto const& base = objectFromHandle(ptr_x); + SELECT_DATA_CASE(base.data_type(), compute_sum_templ, ptr_x, ptr_z); return new DataHandle; } CATCH; } +template +void +compute_max_templ(const void* ptr_x, void* ptr_z) +{ + auto const& x = objectFromHandle >(ptr_x); + *static_cast(ptr_z) = x.max(); +} + extern "C" void* cSIRF_compute_max(const void* ptr_x, void* ptr_z) { try { - auto const& x = objectFromHandle(ptr_x); - x.max(ptr_z); + auto const& base = objectFromHandle(ptr_x); + SELECT_DATA_CASE(base.data_type(), compute_max_templ, ptr_x, ptr_z); return new DataHandle; } CATCH; } +template +void +axpby_templ( + const void* ptr_a, const void* ptr_x, + const void* ptr_b, const void* ptr_y, void* h +) { + auto const& x = objectFromHandle >(ptr_x); + auto const& y = objectFromHandle >(ptr_y); + T a = *static_cast(ptr_a); + T b = *static_cast(ptr_b); + auto& z = objectFromHandle >(h); + z.xapyb(x, a, y, b); +} + extern "C" void* cSIRF_axpby( @@ -185,11 +234,9 @@ cSIRF_axpby( const void* ptr_b, const void* ptr_y ) { try { - auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); + auto const& x = objectFromHandle(ptr_x); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - z.xapyb(x, ptr_a, y, ptr_b); + SELECT_DATA_CASE(x.data_type(), axpby_templ, ptr_a, ptr_x, ptr_b, ptr_y, h); return h; } CATCH; @@ -197,21 +244,33 @@ cSIRF_axpby( extern "C" void* -cSIRF_axpbyAlt( +cSIRF_compute_axpby( const void* ptr_a, const void* ptr_x, const void* ptr_b, const void* ptr_y, void* ptr_z ) { try { - auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); - auto& z = objectFromHandle(ptr_z); - z.xapyb(x, ptr_a, y, ptr_b); + auto const& x = objectFromHandle(ptr_x); + SELECT_DATA_CASE(x.data_type(), axpby_templ, ptr_a, ptr_x, ptr_b, ptr_y, ptr_z); return new DataHandle; } CATCH; } +template +void +xapyb_templ( + const void* ptr_x, const void* ptr_a, + const void* ptr_y, const void* ptr_b, void* h +) { + auto const& x = objectFromHandle >(ptr_x); + auto const& a = objectFromHandle >(ptr_a); + auto const& y = objectFromHandle >(ptr_y); + auto const& b = objectFromHandle >(ptr_b); + auto& z = objectFromHandle >(h); + z.xapyb(x, a, y, b); +} + extern "C" void* cSIRF_xapyb( @@ -220,12 +279,8 @@ cSIRF_xapyb( ) { try { auto const& x = objectFromHandle(ptr_x); - auto const& a = objectFromHandle(ptr_a); - auto const& y = objectFromHandle(ptr_y); - auto const& b = objectFromHandle(ptr_b); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - z.xapyb(x, a, y, b); + SELECT_DATA_CASE(x.data_type(), xapyb_templ, ptr_x, ptr_a, ptr_y, ptr_b, h); return h; } CATCH; @@ -233,23 +288,33 @@ cSIRF_xapyb( extern "C" void* -cSIRF_xapybAlt( +cSIRF_compute_xapyb( const void* ptr_x, const void* ptr_a, const void* ptr_y, const void* ptr_b, void* ptr_z ) { try { - auto const& x = objectFromHandle(ptr_x); - auto const& a = objectFromHandle(ptr_a); - auto const& y = objectFromHandle(ptr_y); - auto const& b = objectFromHandle(ptr_b); - auto& z = objectFromHandle(ptr_z); - z.xapyb(x, a, y, b); + auto const& x = objectFromHandle(ptr_x); + SELECT_DATA_CASE(x.data_type(), xapyb_templ, ptr_x, ptr_a, ptr_y, ptr_b, ptr_z); return new DataHandle; } CATCH; } +template +void +XapYB_templ( + const void* ptr_x, const void* ptr_a, + const void* ptr_y, const void* ptr_b, void* h +) { + auto const& x = objectFromHandle >(ptr_x); + auto const& y = objectFromHandle >(ptr_y); + T a = *static_cast(ptr_a); + auto const& b = objectFromHandle >(ptr_b); + auto& z = objectFromHandle >(h); + z.xapyb(x, a, y, b); +} + extern "C" void* cSIRF_XapYB( @@ -258,11 +323,8 @@ cSIRF_XapYB( ) { try { auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); - auto const& b = objectFromHandle(ptr_b); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - z.xapyb(x, ptr_a, y, b); + SELECT_DATA_CASE(x.data_type(), XapYB_templ, ptr_x, ptr_a, ptr_y, ptr_b, h); return h; } CATCH; @@ -270,30 +332,36 @@ cSIRF_XapYB( extern "C" void* -cSIRF_XapYBAlt( +cSIRF_compute_XapYB( const void* ptr_x, const void* ptr_a, const void* ptr_y, const void* ptr_b, void* ptr_z ) { try { auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); - auto const& b = objectFromHandle(ptr_b); - auto& z = objectFromHandle(ptr_z); - z.xapyb(x, ptr_a, y, b); + SELECT_DATA_CASE(x.data_type(), XapYB_templ, ptr_x, ptr_a, ptr_y, ptr_b, ptr_z); return new DataHandle; } CATCH; } +template +void +add_templ(const void* ptr_x, const void* ptr_y, void* ptr_z) +{ + auto const& x = objectFromHandle >(ptr_x); + T y = *static_cast(ptr_y); + auto& z = objectFromHandle >(ptr_z); + z.add(x, y); +} + extern "C" void* -cSIRF_add(const void* ptr_x, const void* ptr_y, const void* ptr_z) +cSIRF_add(const void* ptr_x, const void* ptr_y, void* ptr_z) { try { auto const& x = objectFromHandle(ptr_x); - auto& z = objectFromHandle(ptr_z); - z.add(x, ptr_y); + SELECT_DATA_CASE(x.data_type(), add_templ, ptr_x, ptr_y, ptr_z); return new DataHandle; } CATCH; @@ -306,34 +374,41 @@ cSIRF_sum(const void* ptr_x, const void* ptr_y) try { auto const& x = objectFromHandle(ptr_x); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - z.add(x, ptr_y); + SELECT_DATA_CASE(x.data_type(), add_templ, ptr_x, ptr_y, h); return h; } CATCH; } +template +void +binary_templ(const void* ptr_x, const void* ptr_y, const char* f, void* h) +{ + auto const& x = objectFromHandle >(ptr_x); + auto const& y = objectFromHandle >(ptr_y); + auto& z = objectFromHandle >(h); + if (sirf::iequals(f, "power")) + z.power(x, y); + else if (sirf::iequals(f, "multiply")) + z.multiply(x, y); + else if (sirf::iequals(f, "divide")) + z.divide(x, y); + else if (sirf::iequals(f, "maximum")) + z.maximum(x, y); + else if (sirf::iequals(f, "minimum")) + z.minimum(x, y); + //else + // return unknownObject("function", f, __FILE__, __LINE__); +} + extern "C" void* cSIRF_binary(const void* ptr_x, const void* ptr_y, const char* f) { try { auto const& x = objectFromHandle(ptr_x); - auto const& y = objectFromHandle(ptr_y); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - if (sirf::iequals(f, "power")) - z.power(x, y); - else if (sirf::iequals(f, "multiply")) - z.multiply(x, y); - else if (sirf::iequals(f, "divide")) - z.divide(x, y); - else if (sirf::iequals(f, "maximum")) - z.maximum(x, y); - else if (sirf::iequals(f, "minimum")) - z.minimum(x, y); - else - return unknownObject("function", f, __FILE__, __LINE__); + SELECT_DATA_CASE(x.data_type(), binary_templ, ptr_x, ptr_y, f, h); return h; } CATCH; @@ -341,29 +416,33 @@ cSIRF_binary(const void* ptr_x, const void* ptr_y, const char* f) extern "C" void* -cSIRF_compute_binary(const void* ptr_x, const void* ptr_y, const char* f, const void* ptr_z) +cSIRF_compute_binary(const void* ptr_x, const void* ptr_y, const char* f, void* ptr_z) { try { auto const& x = objectFromHandle(ptr_x); - auto& z = objectFromHandle(ptr_z); - auto const& y = objectFromHandle(ptr_y); - if (sirf::iequals(f, "power")) - z.power(x, y); - else if (sirf::iequals(f, "multiply")) - z.multiply(x, y); - else if (sirf::iequals(f, "divide")) - z.divide(x, y); - else if (sirf::iequals(f, "maximum")) - z.maximum(x, y); - else if (sirf::iequals(f, "minimum")) - z.minimum(x, y); - else - return unknownObject("function", f, __FILE__, __LINE__); + SELECT_DATA_CASE(x.data_type(), binary_templ, ptr_x, ptr_y, f, ptr_z); return new DataHandle; } CATCH; } +template +void +semibinary_templ(const void* ptr_x, const void* ptr_y, const char* f, void* h) +{ + auto const& x = objectFromHandle >(ptr_x); + auto& z = objectFromHandle >(h); + T y = *static_cast(ptr_y); + if (sirf::iequals(f, "power")) + z.power(x, y); + else if (sirf::iequals(f, "multiply")) + z.multiply(x, y); + else if (sirf::iequals(f, "maximum")) + z.maximum(x, y); + else if (sirf::iequals(f, "minimum")) + z.minimum(x, y); +} + extern "C" void* cSIRF_semibinary(const void* ptr_x, const void* ptr_y, const char* f) @@ -371,17 +450,7 @@ cSIRF_semibinary(const void* ptr_x, const void* ptr_y, const char* f) try { auto const& x = objectFromHandle(ptr_x); void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); - if (sirf::iequals(f, "power")) - z.power(x, ptr_y); - else if (sirf::iequals(f, "multiply")) - z.multiply(x, ptr_y); - else if (sirf::iequals(f, "maximum")) - z.maximum(x, ptr_y); - else if (sirf::iequals(f, "minimum")) - z.minimum(x, ptr_y); - else - return unknownObject("function", f, __FILE__, __LINE__); + SELECT_DATA_CASE(x.data_type(), semibinary_templ, ptr_x, ptr_y, f, h); return h; } CATCH; @@ -389,34 +458,22 @@ cSIRF_semibinary(const void* ptr_x, const void* ptr_y, const char* f) extern "C" void* -cSIRF_compute_semibinary(const void* ptr_x, const void* ptr_y, const char* f, const void* ptr_z) +cSIRF_compute_semibinary(const void* ptr_x, const void* ptr_y, const char* f, void* ptr_z) { try { auto const& x = objectFromHandle(ptr_x); - auto& z = objectFromHandle(ptr_z); - if (sirf::iequals(f, "power")) - z.power(x, ptr_y); - else if (sirf::iequals(f, "multiply")) - z.multiply(x, ptr_y); - else if (sirf::iequals(f, "maximum")) - z.maximum(x, ptr_y); - else if (sirf::iequals(f, "minimum")) - z.minimum(x, ptr_y); - else - return unknownObject("function", f, __FILE__, __LINE__); + SELECT_DATA_CASE(x.data_type(), semibinary_templ, ptr_x, ptr_y, f, ptr_z); return new DataHandle; } CATCH; } -extern "C" -void* -cSIRF_unary(const void* ptr_x, const char* f) +template +void +unary_templ(const void* ptr_x, const char* f, void* h) { - try { - auto const& x = objectFromHandle(ptr_x); - void* h = x.new_data_container_handle(); - auto& z = objectFromHandle(h); + auto const& x = objectFromHandle >(ptr_x); + auto& z = objectFromHandle >(h); if (sirf::iequals(f, "exp")) z.exp(x); else if (sirf::iequals(f, "log")) @@ -427,8 +484,18 @@ cSIRF_unary(const void* ptr_x, const char* f) z.sign(x); else if (sirf::iequals(f, "abs")) z.abs(x); - else - return unknownObject("function", f, __FILE__, __LINE__); + //else + // return unknownObject("function", f, __FILE__, __LINE__); +} + +extern "C" +void* +cSIRF_unary(const void* ptr_x, const char* f) +{ + try { + auto const& x = objectFromHandle(ptr_x); + void* h = x.new_data_container_handle(); + SELECT_DATA_CASE(x.data_type(), unary_templ, ptr_x, f, h); return h; } CATCH; @@ -436,23 +503,11 @@ cSIRF_unary(const void* ptr_x, const char* f) extern "C" void* -cSIRF_compute_unary(const void* ptr_x, const char* f, const void* ptr_z) +cSIRF_compute_unary(const void* ptr_x, const char* f, void* ptr_z) { try { auto const& x = objectFromHandle(ptr_x); - auto& z = objectFromHandle(ptr_z); - if (sirf::iequals(f, "exp")) - z.exp(x); - else if (sirf::iequals(f, "log")) - z.log(x); - else if (sirf::iequals(f, "sqrt")) - z.sqrt(x); - else if (sirf::iequals(f, "sign")) - z.sign(x); - else if (sirf::iequals(f, "abs")) - z.abs(x); - else - return unknownObject("function", f, __FILE__, __LINE__); + SELECT_DATA_CASE(x.data_type(), unary_templ, ptr_x, f, ptr_z); return new DataHandle; } CATCH; @@ -491,14 +546,22 @@ cSIRF_DataHandleVector_push_back(void* self, void* to_append) return new DataHandle; } +template +void +fill_templ(void* ptr_im, const void* ptr_src) +{ + auto& id = objectFromHandle >(ptr_im); + auto const& id_src = objectFromHandle >(ptr_src); + id.fill(id_src); +} + extern "C" void* cSIRF_fillImageFromImage(void* ptr_im, const void* ptr_src) { try { - auto& id = objectFromHandle(ptr_im); - auto const& id_src = objectFromHandle(ptr_src); - id.fill(id_src); + auto& id = objectFromHandle(ptr_im); + SELECT_DATA_CASE(id.data_type(), fill_templ, ptr_im, ptr_src); return new DataHandle; } CATCH; @@ -510,46 +573,74 @@ cSIRF_readImageData(const char* file, const char* eng, int verb) { try { ImageDataWrap idw(file, eng, verb); - std::shared_ptr sptr_id = idw.data_sptr(); - return newObjectHandle(sptr_id); + std::shared_ptr sptr_id = idw.data_sptr(); + return newObjectHandle(sptr_id); } CATCH; } +template +void +equal_images_templ +(const void* ptr_im_a, const void* ptr_im_b, int* ptr_same) +{ + auto const& id_a = objectFromHandle >(ptr_im_a); + auto const& id_b = objectFromHandle >(ptr_im_b); + *ptr_same = (id_a == id_b); +} + extern "C" void* cSIRF_equalImages(const void* ptr_im_a, const void* ptr_im_b) { try { - auto const& id_a = objectFromHandle(ptr_im_a); - auto const& id_b = objectFromHandle(ptr_im_b); - int same = (id_a == id_b); + auto const& id_a = objectFromHandle(ptr_im_a); + int same; + SELECT_DATA_CASE(id_a.data_type(), equal_images_templ, ptr_im_a, ptr_im_b, &same); return dataHandle(same); } CATCH; } +template +void +reorient_templ(void* im_ptr, void* geom_info_ptr) +{ + auto& id = objectFromHandle >(im_ptr); + VoxelisedGeometricalInfo3D geom_info = + objectFromHandle(geom_info_ptr); + id.reorient(geom_info); +} + extern "C" void* cSIRF_ImageData_reorient(void* im_ptr, void *geom_info_ptr) { try { - auto& id = objectFromHandle(im_ptr); - VoxelisedGeometricalInfo3D geom_info = - objectFromHandle(geom_info_ptr); - id.reorient(geom_info); + auto& id = objectFromHandle(im_ptr); + SELECT_DATA_CASE(id.data_type(), reorient_templ, im_ptr, geom_info_ptr); return new DataHandle; } CATCH; } +template +void +geom_info_templ(const void* ptr_im, void** ptr) +{ + const auto& id = objectFromHandle >(ptr_im); + *ptr = newObjectHandle(id.get_geom_info_sptr()); +} + extern "C" void* cSIRF_ImageData_get_geom_info(const void* ptr_im) { try { - const auto& id = objectFromHandle(ptr_im); - return newObjectHandle(id.get_geom_info_sptr()); + const auto& id = objectFromHandle(ptr_im); + void* v; + SELECT_DATA_CASE(id.data_type(), geom_info_templ, ptr_im, &v); + return v; } CATCH; } diff --git a/src/common/include/sirf/Syn/utilities.h b/src/common/include/sirf/Syn/utilities.h index 9e1b0a781..e727ed7ce 100644 --- a/src/common/include/sirf/Syn/utilities.h +++ b/src/common/include/sirf/Syn/utilities.h @@ -34,23 +34,23 @@ namespace sirf { class ImageDataWrap { public: ImageDataWrap(const std::string &filename, const std::string &engine, bool verbose); - ImageData& data() + DataContainer& data() { return *img_sptr_; } - const ImageData& data() const + const DataContainer& data() const { return *img_sptr_; } - std::shared_ptr data_sptr() + std::shared_ptr data_sptr() { return img_sptr_; } - const std::shared_ptr data_sptr() const + const std::shared_ptr data_sptr() const { return img_sptr_; } private: - std::shared_ptr img_sptr_; + std::shared_ptr img_sptr_; }; } diff --git a/src/common/include/sirf/common/DataContainer.h b/src/common/include/sirf/common/DataContainer.h index 67360c164..46124e1de 100644 --- a/src/common/include/sirf/common/DataContainer.h +++ b/src/common/include/sirf/common/DataContainer.h @@ -37,15 +37,22 @@ which rely on the same features of the items. namespace sirf { - typedef std::map Dimensions; - class DataContainer { public: virtual ~DataContainer() {} - //virtual DataContainer* new_data_container() const = 0; + virtual ObjectHandle* new_data_container_handle() const = 0; + + virtual std::string data_type() const = 0; + virtual unsigned int items() const = 0; - virtual bool is_complex() const = 0; + + virtual bool is_complex() const + { + // default value + return false; + } + /// returns the size of data elements virtual int bits() const { @@ -53,73 +60,161 @@ namespace sirf { return is_complex() ? 16 * sizeof(float) : 8 * sizeof(float); } + virtual void write(const std::string& filename) const = 0; + + bool is_empty() const + { + return items() < 1; + } + + std::unique_ptr clone() const + { + return std::unique_ptr(this->clone_impl()); + } + + /// overwrites this container's complex data with complex conjugate values + void conjugate() + { + this->conjugate_impl(); + } + + /// returns unique pointer to the complex-conjugated copy of this container + std::unique_ptr conjugate() const + { + DataContainer* ptr = this->clone_impl(); + ptr->conjugate(); + return std::unique_ptr(ptr); + } + /// returns the norm of this container viewed as a vector virtual float norm() const = 0; - /// below all void* are actually either float* (STIR containers and NiftiImageData) - /// or complex_float_t* (Gadgetron containers) + protected: + virtual DataContainer* clone_impl() const = 0; + /// we assume data to be real, complex data containers must override this + virtual void conjugate_impl() + { + if (is_complex()) + THROW("complex data containes must override conjugate_impl()"); + } + }; - /// calculates the dot product of this container with another one - virtual void dot(const DataContainer& dc, void* ptr) const = 0; + template + class DataContainerTempl : public DataContainer { + public: + virtual ~DataContainerTempl() {} + + virtual ObjectHandle* new_data_container_handle() const = 0; + + virtual std::string data_type() const = 0; + + virtual unsigned int items() const = 0; + + /// returns the norm of this container viewed as a vector + virtual float norm() const = 0; /// calculates the sum of this container elements - virtual void sum(void* ptr) const = 0; + virtual T sum() const = 0; /// calculates the value of this container's element with the largest real part - virtual void max(void* ptr) const = 0; + virtual T max() const = 0; + + /// calculates the dot product of this container with another one + virtual T dot(const DataContainer& dc) const = 0; /// \c *this = the elementwise product \c x*y - virtual void multiply - (const DataContainer& x, const DataContainer& y) = 0; + virtual void multiply(const DataContainer& x, const DataContainer& y) + { + binary_op(x, y, product); + } /// \c *this = the product \c x * y with scalar y - virtual void multiply - (const DataContainer& x, const void* ptr_y) = 0; + void multiply(const DataContainer& x, T y) + { + semibinary_op(x, y, product); + } /// \c *this = the sum \c x + y with scalar y - virtual void add - (const DataContainer& x, const void* ptr_y) = 0; + void add(const DataContainer& x, T y) + { + semibinary_op(x, y, sum); + } /// \c *this = the elementwise ratio \c x / y - virtual void divide - (const DataContainer& x, const DataContainer& y) = 0; + virtual void divide(const DataContainer& x, const DataContainer& y) + { + binary_op(x, y, ratio); + } /// \c *this = the elementwise \c max(x, y) - virtual void maximum - (const DataContainer& x, const DataContainer& y) = 0; - virtual void maximum - (const DataContainer& x, const void* ptr_y) = 0; + void maximum(const DataContainer& x, const DataContainer& y) + { + binary_op(x, y, maxreal); + } + void maximum(const DataContainer& x, T y) + { + semibinary_op(x, y, maxreal); + } /// \c *this = the elementwise \c min(x, y) - virtual void minimum - (const DataContainer& x, const DataContainer& y) = 0; - virtual void minimum - (const DataContainer& x, const void* ptr_y) = 0; + void minimum(const DataContainer& x, const DataContainer& y) + { + binary_op(x, y, minreal); + } + void minimum(const DataContainer& x, T y) + { + semibinary_op(x, y, minreal); + } /// \c *this = the elementwise \c pow(x, y) - virtual void power - (const DataContainer& x, const DataContainer& y) = 0; - virtual void power - (const DataContainer& x, const void* ptr_y) = 0; + void power(const DataContainer& x, const DataContainer& y) + { + binary_op(x, y, power); + //binary_op(x, y, std::pow); + } + void power(const DataContainer& x, T y) + { + semibinary_op(x, y, power); + } /// \c *this = the elementwise \c exp(x) - virtual void exp(const DataContainer& x) = 0; + void exp(const DataContainer& x) + { + unary_op(x, exp); + } /// \c *this = the elementwise \c log(x) - virtual void log(const DataContainer& x) = 0; + void log(const DataContainer& x) + { + unary_op(x, log); + } /// \c *this = the elementwise \c sqrt(x) - virtual void sqrt(const DataContainer& x) = 0; + void sqrt(const DataContainer& x) + { + unary_op(x, sqrt); + } /// \c *this = the elementwise \c sign(x) - virtual void sign(const DataContainer& x) = 0; + void sign(const DataContainer& x) + { + unary_op(x, sign); + } /// \c *this = the elementwise \c abs(x) - virtual void abs(const DataContainer& x) = 0; - + void abs(const DataContainer& x) + { + unary_op(x, abs); + } /// \c *this = the linear combination of \c x and \c y virtual void axpby( - const void* ptr_a, const DataContainer& x, - const void* ptr_b, const DataContainer& y) = 0; - /// alternative interface to the above + T a, const DataContainer& x, + T b, const DataContainer& y) = 0; virtual void xapyb( - const DataContainer& x, const void* ptr_a, - const DataContainer& y, const void* ptr_b) = 0; + const DataContainer& x, T a, + const DataContainer& y, T b) = 0; + //virtual void axpby( + // const void* ptr_a, const DataContainer& x, + // const void* ptr_b, const DataContainer& y) = 0; + /// alternative interface to the above + //virtual void xapyb( + // const DataContainer& x, const void* ptr_a, + // const DataContainer& y, const void* ptr_b) = 0; /// \c *this = elementwise sum of two elementwise products \c x*a and \c y*b virtual void xapyb( @@ -128,134 +223,104 @@ namespace sirf { /// \c *this = elementwise sum of \c x*a and elementwise \c y*b virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, T a, const DataContainer& a_y, const DataContainer& a_b) = 0; + //virtual void xapyb( + // const DataContainer& a_x, const void* ptr_a, + // const DataContainer& a_y, const DataContainer& a_b) = 0; /// \c *this = elementwise sum of elementwise \c x*a and \c y*b void xapyb( const DataContainer& a_x, const DataContainer& a_a, - const DataContainer& a_y, const void* ptr_b) - { - xapyb(a_y, ptr_b, a_x, a_a); - } - - virtual void write(const std::string &filename) const = 0; - - bool is_empty() const + const DataContainer& a_y, T b) { - return items() < 1; + xapyb(a_y, b, a_x, a_a); } + //void xapyb( + // const DataContainer& a_x, const DataContainer& a_a, + // const DataContainer& a_y, const void* ptr_b) + //{ + // xapyb(a_y, ptr_b, a_x, a_a); + //} - std::unique_ptr clone() const - { - return std::unique_ptr(this->clone_impl()); - } + protected: - /// overwrites this container's complex data with complex conjugate values - void conjugate() - { - this->conjugate_impl(); - } + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, T(*f)(T, T)) = 0; + virtual void semibinary_op(const DataContainer& a_x, T y, T(*f)(T, T)) = 0; + virtual void unary_op(const DataContainer& a_x, T(*f)(T)) = 0; - /// returns unique pointer to the complex-conjugated copy of this container - std::unique_ptr conjugate() const - { - DataContainer* ptr = this->clone_impl(); - ptr->conjugate(); - return std::unique_ptr(ptr); - } - - template static T product(T x, T y) { return x * y; } - template static T ratio(T x, T y) { return x / y; } - template static T inverse_ratio(T x, T y) { return y / x; } - template static T sum(T x, T y) { return x + y; } - template static T maximum(T x, T y) { return std::max(x, y); } - template static T maxabs(T x, T y) { return std::max(std::abs(x), std::abs(y)); } - template static T maxreal(T x, T y) { return std::real(x) > std::real(y) ? x : y; } - template static T minimum(T x, T y) { return std::min(x, y); } - template static T minabs(T x, T y) { return std::min(std::abs(x), std::abs(y)); } - template static T minreal(T x, T y) { return std::real(x) < std::real(y) ? x : y; } - static std::complex power(std::complex x, std::complex y) + static T power(T x, T y) { return std::pow(x, y); } - static std::complex exp(std::complex x) + static T exp(T x) { return std::exp(x); } - static std::complex log(std::complex x) + static T log(T x) { return std::log(x); } - static std::complex sqrt(std::complex x) + static T sqrt(T x) { - return std::sqrt(x); + return T(std::sqrt(x)); } - template static T sign(T x) { return (std::real(x) > 0) - (std::real(x) < 0); } - template static T abs(T x) { return T(std::abs(x)); } - - protected: - virtual DataContainer* clone_impl() const = 0; - /// we assume data to be real, complex data containers must override this - virtual void conjugate_impl() - { - if (is_complex()) - THROW("complex data containes must override conjugate_impl()"); - } }; + + typedef std::map Dimensions; } #endif diff --git a/src/common/include/sirf/common/ImageData.h b/src/common/include/sirf/common/ImageData.h index d3e533cd6..9b9383853 100644 --- a/src/common/include/sirf/common/ImageData.h +++ b/src/common/include/sirf/common/ImageData.h @@ -33,8 +33,10 @@ limitations under the License. \brief Abstract base class for SIRF image data. */ + namespace sirf { - class ImageData : public DataContainer + template + class ImageData : public DataContainerTempl { public: virtual ~ImageData() {} @@ -112,7 +114,7 @@ namespace sirf { return !(*this == id); } /// Get geometrical info - std::shared_ptr get_geom_info_sptr() const + virtual std::shared_ptr get_geom_info_sptr() const { // If the geometrical info has not been created yet, throw an error if (!_geom_info_sptr) @@ -128,12 +130,24 @@ namespace sirf { return std::unique_ptr(this->clone_impl()); } /// Is complex? Unless overwridden (Gadgetron), assume not complex. - virtual bool is_complex() const { return false; } + virtual bool is_complex() const { return false; } /// Reorient image. Requires that dimesions and spacing match - virtual void reorient(const VoxelisedGeometricalInfo3D &); - /// Can reorient? (check dimensions and spacing) - static bool can_reorient(const VoxelisedGeometricalInfo3D &geom_1, const VoxelisedGeometricalInfo3D &geom_2, const bool throw_error); - /// Populate the geometrical info metadata (from the image's own metadata) + virtual void reorient(const VoxelisedGeometricalInfo3D &) //; + { + throw std::runtime_error("ImageData::reorient not yet implemented for your image type."); + } + /// Can reorient? (check dimensions and spacing) + static bool can_reorient(const VoxelisedGeometricalInfo3D &geom_1, const VoxelisedGeometricalInfo3D &geom_2, const bool throw_error) //; + { + // If size and spacing match, return true + if (geom_1.get_size() == geom_2.get_size()) + return true; + // Else (and error desired), print error + if (throw_error) + throw std::runtime_error("ImageData::can_reorient: num voxels do not match."); + // Else, return false + return false; + } /// Populate the geometrical info metadata (from the image's own metadata) virtual void set_up_geom_info() = 0; protected: /// Clone helper function. Don't use. diff --git a/src/common/include/sirf/common/csirf.h b/src/common/include/sirf/common/csirf.h index 82b5505db..d744059e7 100644 --- a/src/common/include/sirf/common/csirf.h +++ b/src/common/include/sirf/common/csirf.h @@ -50,30 +50,30 @@ void* cSIRF_compute_sum(const void* ptr_x, PTR_FLOAT ptr_z); void* cSIRF_compute_max(const void* ptr_x, PTR_FLOAT ptr_z); void* cSIRF_axpby(const PTR_FLOAT ptr_a, const void* ptr_x, const PTR_FLOAT ptr_b, const void* ptr_y); -void* cSIRF_axpbyAlt(const PTR_FLOAT ptr_a, const void* ptr_x, +void* cSIRF_compute_axpby(const PTR_FLOAT ptr_a, const void* ptr_x, const PTR_FLOAT ptr_b, const void* ptr_y, void* ptr_z); void* cSIRF_xapyb( const void* ptr_x, const void* ptr_a, const void* ptr_y, const void* ptr_b); -void* cSIRF_xapybAlt( +void* cSIRF_compute_xapyb( const void* ptr_x, const void* ptr_a, const void* ptr_y, const void* ptr_b, void* ptr_z); void* cSIRF_XapYB( const void* ptr_x, const PTR_FLOAT ptr_a, const void* ptr_y, const void* ptr_b); -void* cSIRF_XapYBAlt( +void* cSIRF_compute_XapYB( const void* ptr_x, const PTR_FLOAT ptr_a, const void* ptr_y, const void* ptr_b, void* ptr_z); -void* cSIRF_add(const void* ptr_x, PTR_FLOAT ptr_y, const void* ptr_z); +void* cSIRF_add(const void* ptr_x, PTR_FLOAT ptr_y, void* ptr_z); void* cSIRF_sum(const void* ptr_x, PTR_FLOAT ptr_y); void* cSIRF_binary(const void* ptr_x, const void* ptr_y, const char* f); -void* cSIRF_compute_binary(const void* ptr_x, const void* ptr_y, const char* f, const void* ptr_z); +void* cSIRF_compute_binary(const void* ptr_x, const void* ptr_y, const char* f, void* ptr_z); void* cSIRF_semibinary(const void* ptr_x, PTR_FLOAT ptr_y, const char* f); -void* cSIRF_compute_semibinary(const void* ptr_x, PTR_FLOAT ptr_y, const char* f, const void* ptr_z); +void* cSIRF_compute_semibinary(const void* ptr_x, PTR_FLOAT ptr_y, const char* f, void* ptr_z); void* cSIRF_unary(const void* ptr_x, const char* f); -void* cSIRF_compute_unary(const void* ptr_x, const char* f, const void* ptr_z); +void* cSIRF_compute_unary(const void* ptr_x, const char* f, void* ptr_z); void* cSIRF_write(const void* ptr, const char* filename); void* cSIRF_clone(void* ptr_x); diff --git a/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp b/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp index 8257dc548..ed8eeb7d3 100644 --- a/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp +++ b/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp @@ -266,8 +266,6 @@ void MRAcquisitionData::get_kspace_dimensions(std::vector& dims) const dims.push_back(nc); } - - void MRAcquisitionData::get_data(complex_float_t* z, int a) { @@ -302,7 +300,6 @@ MRAcquisitionData::get_data(complex_float_t* z, int a) void MRAcquisitionData::set_user_floats(float const * const z, int const idx) { - if(idx >= ISMRMRD::ISMRMRD_USER_FLOATS) throw LocalisedException("You try to set the user floats of an index higher than available in the memory of ISMRMRDAcquisition. Pass a smaller idx." , __FILE__, __LINE__); @@ -374,6 +371,7 @@ MRAcquisitionData::xapyb } } +/* void MRAcquisitionData::binary_op (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, @@ -417,106 +415,107 @@ void MRAcquisitionData::multiply (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::binary_op(acq_x, acq_y, DataContainer::product); + MRAcquisitionData::binary_op(acq_x, acq_y, DataContainerTempl::product); } void MRAcquisitionData::multiply (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y) { - MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainer::product); + MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainerTempl::product); } void MRAcquisitionData::add (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y) { - MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainer::sum); + MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainerTempl::sum); } void MRAcquisitionData::divide (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::binary_op(acq_x, acq_y, DataContainer::ratio); + MRAcquisitionData::binary_op(acq_x, acq_y, DataContainerTempl::ratio); } void MRAcquisitionData::maximum (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::binary_op(acq_x, acq_y, DataContainer::maxreal); + MRAcquisitionData::binary_op(acq_x, acq_y, DataContainerTempl::maxreal); } void MRAcquisitionData::maximum (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y) { - MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainer::maxreal); + MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainerTempl::maxreal); } void MRAcquisitionData::minimum (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::binary_op(acq_x, acq_y, DataContainer::minreal); + MRAcquisitionData::binary_op(acq_x, acq_y, DataContainerTempl::minreal); } void MRAcquisitionData::minimum (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y) { - MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainer::minreal); + MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainerTempl::minreal); } void MRAcquisitionData::power (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::binary_op(acq_x, acq_y, DataContainer::power); + MRAcquisitionData::binary_op(acq_x, acq_y, DataContainerTempl::power); } void MRAcquisitionData::power (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y) { - MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainer::power); + MRAcquisitionData::semibinary_op(acq_x, acq_y, y, DataContainerTempl::power); } void MRAcquisitionData::exp (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::unary_op(acq_x, acq_y, DataContainer::exp); + MRAcquisitionData::unary_op(acq_x, acq_y, DataContainerTempl::exp); } void MRAcquisitionData::log (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::unary_op(acq_x, acq_y, DataContainer::log); + MRAcquisitionData::unary_op(acq_x, acq_y, DataContainerTempl::log); } void MRAcquisitionData::sqrt (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::unary_op(acq_x, acq_y, DataContainer::sqrt); + MRAcquisitionData::unary_op(acq_x, acq_y, DataContainerTempl::sqrt); } void MRAcquisitionData::sign (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::unary_op(acq_x, acq_y, DataContainer::sign); + MRAcquisitionData::unary_op(acq_x, acq_y, DataContainerTempl::sign); } void MRAcquisitionData::abs (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y) { - MRAcquisitionData::unary_op(acq_x, acq_y, DataContainer::abs); + MRAcquisitionData::unary_op(acq_x, acq_y, DataContainerTempl::abs); } +*/ complex_float_t MRAcquisitionData::dot @@ -570,8 +569,8 @@ MRAcquisitionData::max(const ISMRMRD::Acquisition& acq_a) return z; } -void -MRAcquisitionData::dot(const DataContainer& dc, void* ptr) const +complex_float_t +MRAcquisitionData::dot(const DataContainer& dc) const { SIRF_DYNAMIC_CAST(const MRAcquisitionData, other, dc); int n = number(); @@ -592,12 +591,11 @@ MRAcquisitionData::dot(const DataContainer& dc, void* ptr) const i++; j++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } -void -MRAcquisitionData::sum(void* ptr) const +complex_float_t +MRAcquisitionData::sum() const { int n = number(); complex_float_t z = 0; @@ -610,12 +608,11 @@ MRAcquisitionData::sum(void* ptr) const z += MRAcquisitionData::sum(a); i++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } -void -MRAcquisitionData::max(void* ptr) const +complex_float_t +MRAcquisitionData::max() const { int n = number(); complex_float_t z = 0; @@ -632,21 +629,18 @@ MRAcquisitionData::max(void* ptr) const z = zi; i++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } void MRAcquisitionData::axpby( - const void* ptr_a, const DataContainer& a_x, - const void* ptr_b, const DataContainer& a_y) + complex_float_t a, const DataContainer& a_x, + complex_float_t b, const DataContainer& a_y) { SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); if (!x.sorted() || !y.sorted()) THROW("a*x + b*y cannot be applied to unsorted x or y"); - complex_float_t a = *static_cast(ptr_a); - complex_float_t b = *static_cast(ptr_b); int nx = x.number(); int ny = y.number(); ISMRMRD::Acquisition ax; @@ -686,27 +680,24 @@ MRAcquisitionData::axpby( void MRAcquisitionData::xapyb( - const DataContainer& a_x, const DataContainer& a_a, + const DataContainer& a_x, complex_float_t a, const DataContainer& a_y, const DataContainer& a_b) { SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, a, a_a); SIRF_DYNAMIC_CAST(const MRAcquisitionData, b, a_b); - if (!x.sorted() || !y.sorted() || !a.sorted() || !b.sorted()) + if (!x.sorted() || !y.sorted() || !b.sorted()) THROW("x*a + y*b cannot be applied to unsorted a, b, x or y"); int nx = x.number(); int ny = y.number(); - int na = a.number(); int nb = b.number(); ISMRMRD::Acquisition ax; ISMRMRD::Acquisition ay; - ISMRMRD::Acquisition aa; ISMRMRD::Acquisition ab; ISMRMRD::Acquisition acq; bool isempty = (number() < 1); - for (int ix = 0, iy = 0, ia = 0, ib = 0, k = 0; - ix < nx && iy < ny && ia < na && ib < nb;) { + for (int ix = 0, iy = 0, ib = 0, k = 0; + ix < nx && iy < ny && ib < nb;) { if (!x.get_acquisition(ix, ax)) { std::cout << ix << " ignored (ax)\n"; ix++; @@ -717,11 +708,6 @@ MRAcquisitionData::xapyb( iy++; continue; } - if (!a.get_acquisition(ia, aa)) { - std::cout << ia << " ignored (aa)\n"; - ia++; - continue; - } if (!b.get_acquisition(ib, ab)) { std::cout << ib << " ignored (ab)\n"; ib++; @@ -734,14 +720,13 @@ MRAcquisitionData::xapyb( continue; } } - MRAcquisitionData::xapyb(ax, aa, ay, ab); + MRAcquisitionData::xapyb(ax, a, ay, ab); if (isempty) append_acquisition(ay); else set_acquisition(k, ay); ix++; iy++; - ia++; ib++; k++; } @@ -751,25 +736,27 @@ MRAcquisitionData::xapyb( void MRAcquisitionData::xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b) { SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); + SIRF_DYNAMIC_CAST(const MRAcquisitionData, a, a_a); SIRF_DYNAMIC_CAST(const MRAcquisitionData, b, a_b); - if (!x.sorted() || !y.sorted() || !b.sorted()) + if (!x.sorted() || !y.sorted() || !a.sorted() || !b.sorted()) THROW("x*a + y*b cannot be applied to unsorted a, b, x or y"); - complex_float_t a = *static_cast(ptr_a); int nx = x.number(); int ny = y.number(); + int na = a.number(); int nb = b.number(); ISMRMRD::Acquisition ax; ISMRMRD::Acquisition ay; + ISMRMRD::Acquisition aa; ISMRMRD::Acquisition ab; ISMRMRD::Acquisition acq; bool isempty = (number() < 1); - for (int ix = 0, iy = 0, ib = 0, k = 0; - ix < nx && iy < ny && ib < nb;) { + for (int ix = 0, iy = 0, ia = 0, ib = 0, k = 0; + ix < nx && iy < ny && ia < na && ib < nb;) { if (!x.get_acquisition(ix, ax)) { std::cout << ix << " ignored (ax)\n"; ix++; @@ -780,6 +767,11 @@ MRAcquisitionData::xapyb( iy++; continue; } + if (!a.get_acquisition(ia, aa)) { + std::cout << ia << " ignored (aa)\n"; + ia++; + continue; + } if (!b.get_acquisition(ib, ab)) { std::cout << ib << " ignored (ab)\n"; ib++; @@ -792,13 +784,14 @@ MRAcquisitionData::xapyb( continue; } } - MRAcquisitionData::xapyb(ax, a, ay, ab); + MRAcquisitionData::xapyb(ax, aa, ay, ab); if (isempty) append_acquisition(ay); else set_acquisition(k, ay); ix++; iy++; + ia++; ib++; k++; } @@ -806,171 +799,62 @@ MRAcquisitionData::xapyb( this->organise_kspace(); } -void -MRAcquisitionData::multiply(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - binary_op(x, y, MRAcquisitionData::multiply); -} - -void -MRAcquisitionData::multiply(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, MRAcquisitionData::multiply); -} - -void -MRAcquisitionData::add(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, MRAcquisitionData::add); -} - -void -MRAcquisitionData::divide(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - binary_op(x, y, MRAcquisitionData::divide); -} - -void -MRAcquisitionData::maximum(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - binary_op(x, y, MRAcquisitionData::maximum); -} - -void -MRAcquisitionData::maximum(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, MRAcquisitionData::maximum); -} - -void -MRAcquisitionData::minimum(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - binary_op(x, y, MRAcquisitionData::minimum); -} - -void -MRAcquisitionData::minimum(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, MRAcquisitionData::minimum); -} - -void -MRAcquisitionData::power(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - binary_op(x, y, MRAcquisitionData::power); -} - -void -MRAcquisitionData::power(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, MRAcquisitionData::power); -} - -void -MRAcquisitionData::exp(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - unary_op(x, MRAcquisitionData::exp); -} - -void -MRAcquisitionData::log(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - unary_op(x, MRAcquisitionData::log); -} - -void -MRAcquisitionData::sqrt(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - unary_op(x, MRAcquisitionData::sqrt); -} - -void -MRAcquisitionData::sign(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - unary_op(x, MRAcquisitionData::sign); -} - -void -MRAcquisitionData::abs(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - unary_op(x, MRAcquisitionData::abs); -} - void MRAcquisitionData::binary_op( const DataContainer& a_x, const DataContainer& a_y, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&)) + complex_float_t(*f)(complex_float_t, complex_float_t)) { - SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); - SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); - if (!x.sorted() || !y.sorted()) - THROW("binary algebraic operations cannot be applied to unsorted data"); + SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); + SIRF_DYNAMIC_CAST(const MRAcquisitionData, y, a_y); + if (!x.sorted() || !y.sorted()) + THROW("binary algebraic operations cannot be applied to unsorted data"); - int nx = x.number(); - int ny = y.number(); - ISMRMRD::Acquisition ax; - ISMRMRD::Acquisition ay; - ISMRMRD::Acquisition acq; - bool isempty = (number() < 1); + int nx = x.number(); + int ny = y.number(); + ISMRMRD::Acquisition ax; + ISMRMRD::Acquisition ay; + ISMRMRD::Acquisition acq; + bool isempty = (number() < 1); for (int ix = 0, iy = 0, k = 0; ix < nx && iy < ny;) { - if (!x.get_acquisition(ix, ax)) { - std::cout << ix << " ignored (ax)\n"; - ix++; - continue; - } - if (!y.get_acquisition(iy, ay)) { - std::cout << iy << " ignored (ay)\n"; - iy++; - continue; - } - if (!isempty) { - if (!get_acquisition(k, acq)) { - std::cout << k << " ignored (acq)\n"; - k++; - continue; - } - } - f(ax, ay); - if (isempty) - append_acquisition(ay); - else - set_acquisition(k, ay); - ix++; - iy++; - k++; - } - this->set_sorted(true); - this->organise_kspace(); + if (!x.get_acquisition(ix, ax)) { + std::cout << ix << " ignored (ax)\n"; + ix++; + continue; + } + if (!y.get_acquisition(iy, ay)) { + std::cout << iy << " ignored (ay)\n"; + iy++; + continue; + } + if (!isempty) { + if (!get_acquisition(k, acq)) { + std::cout << k << " ignored (acq)\n"; + k++; + continue; + } + } + const complex_float_t* px; + complex_float_t* py; + for (px = ax.data_begin(), py = ay.data_begin(); + px != ax.data_end() && py != ay.data_end(); px++, py++) { + *py = f(*px, *py); + } + if (isempty) + append_acquisition(ay); + else + set_acquisition(k, ay); + ix++; + iy++; + k++; + } + this->set_sorted(true); + this->organise_kspace(); } void -MRAcquisitionData::semibinary_op(const DataContainer& a_x, complex_float_t y, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&, complex_float_t y)) +MRAcquisitionData::semibinary_op( + const DataContainer& a_x, complex_float_t y, + complex_float_t(*f)(complex_float_t, complex_float_t)) { SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); if (!x.sorted()) @@ -995,7 +879,12 @@ MRAcquisitionData::semibinary_op(const DataContainer& a_x, complex_float_t y, } } x.get_acquisition(ix, ay); - f(ax, ay, y); + const complex_float_t* px; + complex_float_t* py; + for (px = ax.data_begin(), py = ay.data_begin(); + px != ax.data_end() && py != ay.data_end(); px++, py++) { + *py = f(*px, y); + } if (isempty) append_acquisition(ay); else @@ -1008,8 +897,9 @@ MRAcquisitionData::semibinary_op(const DataContainer& a_x, complex_float_t y, } void -MRAcquisitionData::unary_op(const DataContainer& a_x, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&)) +MRAcquisitionData::unary_op( + const DataContainer& a_x, + complex_float_t(*f)(complex_float_t)) { SIRF_DYNAMIC_CAST(const MRAcquisitionData, x, a_x); if (!x.sorted()) @@ -1034,7 +924,12 @@ MRAcquisitionData::unary_op(const DataContainer& a_x, } } x.get_acquisition(ix, ay); - f(ax, ay); + const complex_float_t* px; + complex_float_t* py; + for (px = ax.data_begin(), py = ay.data_begin(); + px != ax.data_end() && py != ay.data_end(); px++, py++) { + *py = f(*px); + } if (isempty) append_acquisition(ay); else @@ -1062,7 +957,6 @@ MRAcquisitionData::norm() const return std::sqrt(r); } - ISMRMRD::TrajectoryType MRAcquisitionData::get_trajectory_type() const { @@ -1073,7 +967,6 @@ MRAcquisitionData::get_trajectory_type() const std::cout << "You have a file with " << hdr.encoding.size() << " encodings. Just the first one is picked." << std::endl; return hdr.encoding[0].trajectory; - } void MRAcquisitionData::set_trajectory_type(const ISMRMRD::TrajectoryType type) @@ -1130,7 +1023,6 @@ MRAcquisitionData::sort_by_time() this->organise_kspace(); sorted_ = true; - } std::vector MRAcquisitionData::get_kspace_order() const @@ -1139,7 +1031,6 @@ std::vector MRAcquisitionData::get_kspace_order() const throw LocalisedException("Your acquisition data object contains no data, so no order is determined." , __FILE__, __LINE__); else if(this->sorting_.size() == 0) throw LocalisedException("The kspace is not sorted yet. Please call organise_kspace(), sort() or sort_by_time() first." , __FILE__, __LINE__); - std::vector output; for(unsigned i = 0; i MRAcquisitionData::get_slice_encoding_index(const unsigned kspa return slice_encode_index; } - void MRAcquisitionData::get_subset(MRAcquisitionData& subset, const std::vector subset_idx) const { subset.set_acquisitions_info(this->acquisitions_info()); @@ -1382,7 +1272,6 @@ KSpaceSubset::TagType KSpaceSubset::get_tag_from_img(const CFImage& img) return tag; } - KSpaceSubset::TagType KSpaceSubset::get_tag_from_acquisition(ISMRMRD::Acquisition acq) { TagType tag; @@ -1416,8 +1305,8 @@ void KSpaceSubset::print_acquisition_tag(ISMRMRD::Acquisition acq) print_tag(tag); } -void -GadgetronImageData::dot(const DataContainer& dc, void* ptr) const +complex_float_t +GadgetronImageData::dot(const DataContainer& dc) const { SIRF_DYNAMIC_CAST(const GadgetronImageData, ic, dc); complex_float_t z = 0; @@ -1426,12 +1315,11 @@ GadgetronImageData::dot(const DataContainer& dc, void* ptr) const const ImageWrap& v = ic.image_wrap(i); z += u.dot(v); } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } -void -GadgetronImageData::sum(void* ptr) const +complex_float_t +GadgetronImageData::sum() const { complex_float_t z = 0; for (unsigned int i = 0; i < number(); i++) { @@ -1439,12 +1327,11 @@ GadgetronImageData::sum(void* ptr) const complex_float_t t = u.sum(); z += t; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } -void -GadgetronImageData::max(void* ptr) const +complex_float_t +GadgetronImageData::max() const { complex_float_t z = 0; for (unsigned int i = 0; i < number(); i++) { @@ -1455,17 +1342,14 @@ GadgetronImageData::max(void* ptr) const if (ri > r) z = zi; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; } void GadgetronImageData::axpby( -const void* ptr_a, const DataContainer& a_x, -const void* ptr_b, const DataContainer& a_y) + complex_float_t a, const DataContainer& a_x, + complex_float_t b, const DataContainer& a_y) { - complex_float_t a = *static_cast(ptr_a); - complex_float_t b = *static_cast(ptr_b); SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); unsigned int nx = x.number(); @@ -1566,127 +1450,6 @@ GadgetronImageData::unary_op(const DataContainer& a_x, this->set_meta_data(x.get_meta_data()); } -void -GadgetronImageData::multiply(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::product); -} - -void -GadgetronImageData::multiply(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::product); -} - -void -GadgetronImageData::add(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::sum); -} - -void -GadgetronImageData::divide(const DataContainer& a_x, const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::ratio); -} - -void -GadgetronImageData::maximum( - const DataContainer& a_x, - const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::maxreal); -} - -void -GadgetronImageData::maximum(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::maxreal); -} - -void -GadgetronImageData::minimum( - const DataContainer& a_x, - const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::minreal); -} - -void -GadgetronImageData::minimum(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::minreal); -} - -void -GadgetronImageData::power( - const DataContainer& a_x, - const DataContainer& a_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::power); -} - -void -GadgetronImageData::power(const DataContainer& a_x, const void* ptr_y) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - complex_float_t y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::power); -} - -void -GadgetronImageData::exp(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::exp); -} - -void -GadgetronImageData::log(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::log); -} - -void -GadgetronImageData::sqrt(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::sqrt); -} - -void -GadgetronImageData::sign(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::sign); -} - -void -GadgetronImageData::abs(const DataContainer& a_x) -{ - SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::abs); -} - float GadgetronImageData::norm() const { @@ -2025,7 +1788,6 @@ GadgetronImagesVector::GadgetronImagesVector(const MRAcquisitionData& ad, const set_meta_data(ad.acquisitions_info()); } - GadgetronImagesVector::GadgetronImagesVector (const GadgetronImagesVector& images) : images_() @@ -2240,7 +2002,7 @@ GadgetronImagesVector::print_header(const unsigned im_num) std::cout << acqs_info_.c_str() << "\n"; } } - +/* bool GadgetronImagesVector::is_complex() const { // If any of the wraps are complex, return true. for (unsigned i=0; iget_geom_info_sptr(); @@ -2524,7 +2286,6 @@ CoilImagesVector::calculate(const MRAcquisitionData& ad) } } - std::unique_ptr CoilImagesVector::extract_calibration_data(const MRAcquisitionData& ad) const { using ISMRMRD::ISMRMRD_AcquisitionFlags; @@ -2697,7 +2458,6 @@ void CoilSensitivitiesVector::combine_images_with_coilmaps(GadgetronImageData& c } } - void CoilSensitivitiesVector::calculate(CoilImagesVector& iv) { @@ -2918,4 +2678,3 @@ CoilSensitivitiesVector::max_diff_ } return s; } - diff --git a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h index 80e30bc66..1c6ceaf31 100644 --- a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h +++ b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h @@ -214,8 +214,15 @@ namespace sirf { \brief Abstract MR acquisition data container class. */ - class MRAcquisitionData : public DataContainer { + class MRAcquisitionData : public DataContainerTempl { public: + typedef MRAcquisitionData self_type; + + virtual std::string data_type() const + { + return std::string("complex_float"); + } + // static methods // ISMRMRD acquisitions algebra: acquisitions viewed as vectors of @@ -250,6 +257,7 @@ namespace sirf { static complex_float_t sum(const ISMRMRD::Acquisition& acq_x); // the value of the element of x with the largest real part static complex_float_t max(const ISMRMRD::Acquisition& acq_x); +/* // elementwise multiplication // y := x .* y static void multiply @@ -299,6 +307,7 @@ namespace sirf { // y := abs(x) static void abs (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y); +*/ static float norm(const ISMRMRD::Acquisition& acq_x); // type and dimension of an ISMRMRD::Acquisition parameter @@ -504,7 +513,7 @@ namespace sirf { return std::make_tuple(limit.minimum, limit.maximum, limit.center); } - // abstract methods + // virtual methods virtual void empty() = 0; virtual void take_over(MRAcquisitionData&) = 0; @@ -537,58 +546,36 @@ namespace sirf { } // acquisition data algebra - /// below all void* are actually complex_float_t* - virtual void sum(void* ptr) const; - virtual void max(void* ptr) const; - virtual void dot(const DataContainer& dc, void* ptr) const; - complex_float_t dot(const DataContainer& a_x) - { - complex_float_t z; - dot(a_x, &z); - return z; - } + virtual float norm() const; + virtual complex_float_t sum() const; + virtual complex_float_t max() const; + virtual complex_float_t dot(const DataContainer& dc) const; virtual void axpby( - const void* ptr_a, const DataContainer& a_x, - const void* ptr_b, const DataContainer& a_y); - virtual void xapyb( - const DataContainer& a_x, const DataContainer& a_a, - const DataContainer& a_y, const DataContainer& a_b); + complex_float_t a, const DataContainer& a_x, + complex_float_t b, const DataContainer& a_y); virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, - const DataContainer& a_y, const void* ptr_b) + const DataContainer& a_x, complex_float_t a, + const DataContainer& a_y, complex_float_t b) { - axpby(ptr_a, a_x, ptr_b, a_y); + axpby(a, a_x, b, a_y); } virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, complex_float_t a, + const DataContainer& a_y, const DataContainer& a_b); + virtual void xapyb( + const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b); - virtual void multiply(const DataContainer& x, const DataContainer& y); - virtual void divide(const DataContainer& x, const DataContainer& y); - virtual void maximum(const DataContainer& x, const DataContainer& y); - virtual void minimum(const DataContainer& x, const DataContainer& y); - virtual void power(const DataContainer& x, const DataContainer& y); - virtual void multiply(const DataContainer& x, const void* y); - virtual void add(const DataContainer& x, const void* ptr_y); - virtual void maximum(const DataContainer& x, const void* y); - virtual void minimum(const DataContainer& x, const void* y); - virtual void power(const DataContainer& x, const void* y); - virtual void exp(const DataContainer& x); - virtual void log(const DataContainer& x); - virtual void sqrt(const DataContainer& x); - virtual void sign(const DataContainer& x); - virtual void abs(const DataContainer& x); - virtual float norm() const; virtual void write(const std::string &filename) const; - // regular methods - void binary_op(const DataContainer& a_x, const DataContainer& a_y, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&)); - void semibinary_op(const DataContainer& a_x, complex_float_t y, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&, complex_float_t)); - void unary_op(const DataContainer& a_x, - void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&)); + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, + complex_float_t(*f)(complex_float_t, complex_float_t)); + virtual void semibinary_op( + const DataContainer& a_x, complex_float_t y, + complex_float_t(*f)(complex_float_t, complex_float_t)); + virtual void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t)); + // regular methods AcquisitionsInfo acquisitions_info() const { return acqs_info_; } void set_acquisitions_info(std::string info) { acqs_info_ = info; } void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;} @@ -716,7 +703,7 @@ namespace sirf { virtual void empty(); virtual void take_over(MRAcquisitionData& ad) {} virtual unsigned int number() const { return (unsigned int)acqs_.size(); } - virtual unsigned int items() const { return (unsigned int)acqs_.size(); } + virtual unsigned int items() const { return (unsigned int)this->acqs_.size(); } virtual void append_acquisition(ISMRMRD::Acquisition& acq) { acqs_.push_back(gadgetron::shared_ptr @@ -782,11 +769,17 @@ namespace sirf { */ - class ISMRMRDImageData : public ImageData { + class ISMRMRDImageData : public ImageData { public: + typedef ISMRMRDImageData self_type; + //ISMRMRDImageData(ISMRMRDImageData& id, const char* attr, //const char* target); //does not build, have to be in the derived class + virtual bool is_complex() const + { + return true; + } virtual void empty() = 0; virtual unsigned int number() const = 0; virtual gadgetron::shared_ptr sptr_image_wrap @@ -885,37 +878,28 @@ namespace sirf { } virtual float norm() const; - /// below all void* are actually complex_float_t* - virtual void sum(void* ptr) const; - virtual void max(void* ptr) const; - virtual void dot(const DataContainer& dc, void* ptr) const; - virtual void axpby( - const void* ptr_a, const DataContainer& a_x, - const void* ptr_b, const DataContainer& a_y); - virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, - const DataContainer& a_y, const void* ptr_b) + virtual complex_float_t sum() const; + virtual complex_float_t max() const; + virtual complex_float_t dot(const DataContainer& dc) const; + void axpby( + complex_float_t a, const DataContainer& a_x, + complex_float_t b, const DataContainer& a_y); + void xapyb( + const DataContainer& a_x, complex_float_t a_a, + const DataContainer& a_y, complex_float_t a_b) { - ComplexFloat_ a(*static_cast(ptr_a)); - ComplexFloat_ b(*static_cast(ptr_b)); + ComplexFloat_ a(a_a); + ComplexFloat_ b(a_b); xapyb_(a_x, a, a_y, b); } virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, complex_float_t a_a, const DataContainer& a_y, const DataContainer& a_b) { - ComplexFloat_ a(*static_cast(ptr_a)); + ComplexFloat_ a(a_a); SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b); xapyb_(a_x, a, a_y, b); } - //virtual void xapyb( - // const DataContainer& a_x, const DataContainer& a_a, - // const DataContainer& a_y, const void* ptr_b) - //{ - // SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a); - // ComplexFloat_ b(*(complex_float_t*)ptr_b); - // xapyb_(a_x, a, a_y, b); - //} virtual void xapyb( const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b) @@ -924,50 +908,17 @@ namespace sirf { SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b); xapyb_(a_x, a, a_y, b); } - virtual void multiply(const DataContainer& x, const DataContainer& y); - virtual void divide(const DataContainer& x, const DataContainer& y); - virtual void maximum(const DataContainer& x, const DataContainer& y); - virtual void minimum(const DataContainer& x, const DataContainer& y); - virtual void power(const DataContainer& x, const DataContainer& y); - virtual void multiply(const DataContainer& x, const void* ptr_y); - virtual void add(const DataContainer& x, const void* ptr_y); - virtual void maximum(const DataContainer& x, const void* ptr_y); - virtual void minimum(const DataContainer& x, const void* ptr_y); - virtual void power(const DataContainer& x, const void* ptr_y); - virtual void exp(const DataContainer& x); - virtual void log(const DataContainer& x); - virtual void sqrt(const DataContainer& x); - virtual void sign(const DataContainer& x); - virtual void abs(const DataContainer& x); - - void binary_op( - const DataContainer& a_x, const DataContainer& a_y, + + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, complex_float_t(*f)(complex_float_t, complex_float_t)); - void semibinary_op( + virtual void semibinary_op( const DataContainer& a_x, complex_float_t y, complex_float_t(*f)(complex_float_t, complex_float_t)); - void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t)); + virtual void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t)); void fill(float s); void scale(float s); - complex_float_t dot(const DataContainer& a_x) - { - complex_float_t z; - dot(a_x, &z); - return z; - } - void axpby( - complex_float_t a, const DataContainer& a_x, - complex_float_t b, const DataContainer& a_y) - { - axpby(&a, a_x, &b, a_y); - } - void xapyb( - const DataContainer& a_x, complex_float_t a, - const DataContainer& a_y, complex_float_t b) - { - xapyb(a_x, &a, a_y, &b); - } + gadgetron::unique_ptr clone() const { return gadgetron::unique_ptr(this->clone_impl()); @@ -1083,6 +1034,10 @@ namespace sirf { class GadgetronImagesVector : public GadgetronImageData { public: + virtual std::string data_type() const + { + return std::string("complex_float"); + } typedef ImageData::Iterator BaseIter; typedef ImageData::Iterator_const BaseIter_const; typedef std::vector >::iterator @@ -1262,7 +1217,7 @@ namespace sirf { } virtual unsigned int items() const { - return (unsigned int)images_.size(); + return (unsigned int)this->images_.size(); } virtual unsigned int number() const { @@ -1385,7 +1340,7 @@ namespace sirf { void print_header(const unsigned im_num); /// Is complex? - virtual bool is_complex() const; + //virtual bool is_complex() const; /// Reorient image. Requires that dimensions match virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out); diff --git a/src/xGadgetron/cGadgetron/tests/mrtests.cpp b/src/xGadgetron/cGadgetron/tests/mrtests.cpp index 1843f20a2..5c92d48bd 100644 --- a/src/xGadgetron/cGadgetron/tests/mrtests.cpp +++ b/src/xGadgetron/cGadgetron/tests/mrtests.cpp @@ -393,13 +393,13 @@ bool test_acq_mod_adjointness(MRAcquisitionData& ad) shared_ptr sptr_random = std::move(sptr_bwd->clone()); sptr_random->set_data(&random_data[0]); - complex_float_t Eh_kdat_Dot_img; - sptr_bwd->dot(*sptr_random, &Eh_kdat_Dot_img); + complex_float_t Eh_kdat_Dot_img = sptr_bwd->dot(*sptr_random); + //sptr_bwd->dot(*sptr_random, &Eh_kdat_Dot_img); auto sptr_fwd = AM.fwd(*sptr_random); - complex_float_t E_img_Dot_kdat; - ad.dot(*sptr_fwd, &E_img_Dot_kdat); + complex_float_t E_img_Dot_kdat = ad.dot(*sptr_fwd); + //ad.dot(*sptr_fwd, &E_img_Dot_kdat); std::cout << "Backward kdata dot random image: " << Eh_kdat_Dot_img << std::endl; std::cout << "Forward random image dot kdata : " << E_img_Dot_kdat << std::endl; diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index 801c78010..b039a875e 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -1241,7 +1241,7 @@ extern "C" void* cSTIR_imageFromImageData(void* ptr_v) { try { - ImageData& id = objectFromHandle(ptr_v); + ImageData& id = objectFromHandle >(ptr_v); shared_ptr sptr(new STIRImageData(id)); return (void*)newObjectHandle(sptr); } diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h index a1499385c..6633f6350 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h @@ -145,9 +145,13 @@ namespace sirf { storage mode (file/memory) selection. */ - class STIRAcquisitionData : public DataContainer { + class STIRAcquisitionData : public DataContainerTempl { public: virtual ~STIRAcquisitionData() {} + virtual std::string data_type() const + { + return std::string("float"); + } // virtual constructors virtual STIRAcquisitionData* same_acquisition_data @@ -246,111 +250,34 @@ namespace sirf { // data container methods unsigned int items() const { - if (_is_empty != -1) - return _is_empty ? 0 : 1; + if (this->_is_empty != -1) + return this->_is_empty ? 0 : 1; try { - get_segment_by_sinogram(0); + this->get_segment_by_sinogram(0); } catch (...) { - _is_empty = 1; + this->_is_empty = 1; return 0; // no data found - this must be a template } - _is_empty = 0; + this->_is_empty = 0; return 1; // data ok } virtual float norm() const; - /// below all void* are actually float* - virtual void sum(void* ptr) const; - virtual void max(void* ptr) const; - virtual void dot(const DataContainer& a_x, void* ptr) const; - float dot(const DataContainer& a_x) const - { - float s; - dot(a_x, &s); - return s; - } + virtual float sum() const; + virtual float max() const; + virtual float dot(const DataContainer& a_x) const; virtual void axpby( - const void* ptr_a, const DataContainer& a_x, - const void* ptr_b, const DataContainer& a_y); + float a, const DataContainer& a_x, + float b, const DataContainer& a_y); virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, - const DataContainer& a_y, const void* ptr_b); + const DataContainer& a_x, float a, + const DataContainer& a_y, float b); virtual void xapyb( - const DataContainer& a_x, const DataContainer& a_a, + const DataContainer& a_x, float a, const DataContainer& a_y, const DataContainer& a_b); virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b); - virtual void abs(const DataContainer& x) - { - unary_op(x, std::abs); - } - virtual void exp(const DataContainer& x) - { - unary_op(x, std::exp); - } - virtual void log(const DataContainer& x) - { - unary_op(x, std::log); - } - virtual void sqrt(const DataContainer& x) - { - unary_op(x, std::sqrt); - } - virtual void sign(const DataContainer& x) - { - unary_op(x, DataContainer::sign); - } - virtual void multiply(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::product); - } - virtual void add(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::sum); - } - virtual void divide(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::ratio); - } - virtual void maximum(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::maximum); - } - virtual void minimum(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::minimum); - } - virtual void power(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, std::pow); - } - virtual void multiply(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::product); - } - virtual void divide(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::ratio); - } - virtual void maximum(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::maximum); - } - virtual void minimum(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::minimum); - } - virtual void power(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, std::pow); - } virtual void inv(float a, const DataContainer& x); virtual void write(const std::string &filename) const { @@ -451,9 +378,11 @@ namespace sirf { (sptr_s, span, max_ring_diff, num_views, num_tang_pos, false)); } - void unary_op(const DataContainer& a_x, float(*f)(float)); - void semibinary_op(const DataContainer& a_x, float y, float(*f)(float, float)); - void binary_op(const DataContainer& a_x, const DataContainer& a_y, float(*f)(float, float)); + //void unary_op(const DataContainer& a_x, float(*f)(float)); + virtual void unary_op(const DataContainer& a_x, float(*f)(float)); + virtual void semibinary_op(const DataContainer& a_x, float y, float(*f)(float, float)); + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, float(*f)(float, float)); + //virtual void binary_op_new(const DataContainer& a_x, const DataContainer& a_y, float(*f)(float, float)) {} protected: static std::string _storage_scheme; @@ -761,7 +690,7 @@ namespace sirf { t += (*iter) * (*iter); return std::sqrt((float)t); } - virtual void dot(const DataContainer& a_x, void* ptr) const + virtual float dot(const DataContainer& a_x) const { auto x = dynamic_cast(&a_x); // Can only do this if both are STIRAcquisitionDataInMemory @@ -769,7 +698,7 @@ namespace sirf { const stir::ProjDataInMemory *pd2_ptr = dynamic_cast(x->data().get()); // If either cast failed, fall back to general method if (is_null_ptr(pd_ptr) || is_null_ptr(pd2_ptr)) - return this->STIRAcquisitionData::dot(a_x,ptr); + return this->STIRAcquisitionData::dot(a_x); // do it double t = 0.0; @@ -777,9 +706,9 @@ namespace sirf { auto iter_other = pd2_ptr->begin(); while (iter != pd_ptr->end()) t += (*iter++) * double(*iter_other++); - - float* ptr_t = static_cast(ptr); - *ptr_t = (float)t; + return (float)t; + //float* ptr_t = static_cast(ptr); + //*ptr_t = (float)t; } virtual void multiply(const DataContainer& x, const DataContainer& y) { @@ -793,7 +722,8 @@ namespace sirf { // If either cast failed, fall back to general method if (is_null_ptr(pd_ptr) || is_null_ptr(pd_x_ptr) || is_null_ptr(pd_x_ptr)) - return this->STIRAcquisitionData::multiply(x,y); + //return this->STIRAcquisitionData::multiply(x, y); + return this->DataContainerTempl::multiply(x, y); // do it auto iter = pd_ptr->begin(); @@ -814,7 +744,8 @@ namespace sirf { // If either cast failed, fall back to general method if (is_null_ptr(pd_ptr) || is_null_ptr(pd_x_ptr) || is_null_ptr(pd_x_ptr)) - return this->STIRAcquisitionData::divide(x,y); + //return this->STIRAcquisitionData::divide(x,y); + return this->DataContainerTempl::divide(x, y); // do it auto iter = pd_ptr->begin(); @@ -852,8 +783,12 @@ namespace sirf { additioanally, implements the linear algebra functionality specified by the abstract base class aDatacontainer. */ - class STIRImageData : public ImageData { + class STIRImageData : public ImageData { public: + virtual std::string data_type() const + { + return std::string("float"); + } typedef ImageData::Iterator BaseIter; typedef ImageData::Iterator_const BaseIter_const; class Iterator : public BaseIter { @@ -1048,92 +983,21 @@ namespace sirf { virtual void write(const std::string& filename, const std::string& format_file) const; virtual float norm() const; - /// below all void* are actually float* - virtual void sum(void* ptr) const; - virtual void max(void* ptr) const; - virtual void dot(const DataContainer& a_x, void* ptr) const; + virtual float sum() const; + virtual float max() const; + virtual float dot(const DataContainer& dc) const; virtual void axpby( - const void* ptr_a, const DataContainer& a_x, - const void* ptr_b, const DataContainer& a_y); + float a, const DataContainer& a_x, + float b, const DataContainer& a_y); virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, - const DataContainer& a_y, const void* ptr_b); + const DataContainer& a_x, float a, + const DataContainer& a_y, float b); virtual void xapyb( - const DataContainer& a_x, const DataContainer& a_a, + const DataContainer& a_x, float a, const DataContainer& a_y, const DataContainer& a_b); virtual void xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, const DataContainer& a_a, const DataContainer& a_y, const DataContainer& a_b); - virtual void abs(const DataContainer& x) - { - unary_op(x, std::abs); - } - virtual void exp(const DataContainer& x) - { - unary_op(x, std::exp); - } - virtual void log(const DataContainer& x) - { - unary_op(x, std::log); - } - virtual void sqrt(const DataContainer& x) - { - unary_op(x, std::sqrt); - } - virtual void sign(const DataContainer& x) - { - unary_op(x, DataContainer::sign); - } - virtual void multiply(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::product); - } - virtual void add(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::sum); - } - virtual void divide(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::ratio); - } - virtual void maximum(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::maximum); - } - virtual void minimum(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::minimum); - } - virtual void power(const DataContainer& x, const void* ptr_y) - { - float y = *static_cast(ptr_y); - semibinary_op(x, y, std::pow); - } - virtual void multiply(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::product); - } - virtual void divide(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::ratio); - } - virtual void maximum(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::maximum); - } - virtual void minimum(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, DataContainer::minimum); - } - virtual void power(const DataContainer& x, const DataContainer& y) - { - binary_op(x, y, std::pow); - } Image3DF& data() { @@ -1169,24 +1033,18 @@ namespace sirf { _data->fill(v); } void scale(float s); - float dot(const DataContainer& a_x) const - { - float s; - dot(a_x, &s); - return s; - } - void axpby( - float a, const DataContainer& a_x, - float b, const DataContainer& a_y) - { - axpby(&a, a_x, &b, a_y); - } - void xapyb( - const DataContainer& a_x, float a, - const DataContainer& a_y, float b) - { - xapyb(a_x, &a, a_y, &b); - } + //void axpby( + // float a, const DataContainer& a_x, + // float b, const DataContainer& a_y) + //{ + // axpby(&a, a_x, &b, a_y); + //} + //void xapyb( + // const DataContainer& a_x, float a, + // const DataContainer& a_y, float b) + //{ + // xapyb(a_x, &a, a_y, &b); + //} size_t size() const { return _data->size_all(); @@ -1252,9 +1110,9 @@ namespace sirf { /// Populate the geometrical info metadata (from the image's own metadata) virtual void set_up_geom_info(); - void unary_op(const DataContainer& a_x, float(*f)(float)); - void semibinary_op(const DataContainer& a_x, float y, float(*f)(float, float)); - void binary_op(const DataContainer& a_x, const DataContainer& a_y, float(*f)(float, float)); + virtual void unary_op(const DataContainer& a_x, float(*f)(float)); + virtual void semibinary_op(const DataContainer& a_x, float y, float(*f)(float, float)); + virtual void binary_op(const DataContainer& a_x, const DataContainer& a_y, float(*f)(float, float)); private: /// Clone helper function. Don't use. diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index 73f7d8f4f..a4d531b24 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -1045,7 +1045,8 @@ The actual algorithm is described in if (have_a) { auto sptr_a = am.additive_term_sptr(); float a = 1.0f; - sptr->axpby(&a, *sptr, &a, *sptr_a); + sptr->axpby(a, *sptr, a, *sptr_a); + //sptr->axpby(&a, *sptr, &a, *sptr_a); } set_additive_proj_data_sptr(sptr->data()); } diff --git a/src/xSTIR/cSTIR/stir_data_containers.cpp b/src/xSTIR/cSTIR/stir_data_containers.cpp index 5211a10c6..590522e1e 100644 --- a/src/xSTIR/cSTIR/stir_data_containers.cpp +++ b/src/xSTIR/cSTIR/stir_data_containers.cpp @@ -55,8 +55,8 @@ STIRAcquisitionData::norm() const return std::sqrt((float)t); } -void -STIRAcquisitionData::sum(void* ptr) const +float +STIRAcquisitionData::sum() const { int n = get_max_segment_num(); double t = 0; @@ -72,12 +72,13 @@ STIRAcquisitionData::sum(void* ptr) const t += *seg_iter++; } } - float* ptr_t = static_cast(ptr); - *ptr_t = (float)t; + return (float)t; + //float* ptr_t = static_cast(ptr); + //*ptr_t = (float)t; } -void -STIRAcquisitionData::max(void* ptr) const +float +STIRAcquisitionData::max() const { int n = get_max_segment_num(); float t = 0; @@ -93,12 +94,13 @@ STIRAcquisitionData::max(void* ptr) const t = std::max(t, *seg_iter++); } } - float* ptr_t = static_cast(ptr); - *ptr_t = (float)t; + return (float)t; + //float* ptr_t = static_cast(ptr); + //*ptr_t = (float)t; } -void -STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const +float +STIRAcquisitionData::dot(const DataContainer& a_x) const { SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x); int n = get_max_segment_num(); @@ -124,39 +126,38 @@ STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const t += (*seg_iter++) * double(*sx_iter++); } } - float* ptr_t = static_cast(ptr); - *ptr_t = (float)t; + return (float)t; + //float* ptr_t = static_cast(ptr); + //*ptr_t = (float)t; } void STIRAcquisitionData::axpby( -const void* ptr_a, const DataContainer& a_x, -const void* ptr_b, const DataContainer& a_y + float a, const DataContainer& a_x, + float b, const DataContainer& a_y ) { //Add deprecation warning - STIRAcquisitionData::xapyb(a_x, ptr_a, a_y, ptr_b); + STIRAcquisitionData::xapyb(a_x, a, a_y, b); } void STIRAcquisitionData::xapyb( -const DataContainer& a_x, const void* ptr_a, -const DataContainer& a_y, const void* ptr_b + const DataContainer& a_x, float a, + const DataContainer& a_y, float b ) { - // Cast to correct types - float a = *static_cast(ptr_a); - float b = *static_cast(ptr_b); + // Cast to correct types auto x = dynamic_cast(&a_x); - auto y = dynamic_cast(&a_y); + auto y = dynamic_cast(&a_y); - if (is_null_ptr(x) || is_null_ptr(x->data()) || - is_null_ptr(y) || is_null_ptr(y->data())) - throw std::runtime_error("STIRAcquisitionData::xapyb: At least one argument is not" - "STIRAcquisitionData or is not initialised."); + if (is_null_ptr(x) || is_null_ptr(x->data()) || + is_null_ptr(y) || is_null_ptr(y->data())) + throw std::runtime_error("STIRAcquisitionData::xapyb: At least one argument is not" + "STIRAcquisitionData or is not initialised."); - // Call STIR's xapyb - data()->xapyb(*x->data(), a, *y->data(), b); + // Call STIR's xapyb + data()->xapyb(*x->data(), a, *y->data(), b); } void @@ -323,13 +324,12 @@ STIRAcquisitionData::binary_op( void STIRAcquisitionData::xapyb( - const DataContainer& a_x, const void* ptr_a, + const DataContainer& a_x, float a, const DataContainer& a_y, const DataContainer& a_b) { SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x); SIRF_DYNAMIC_CAST(const STIRAcquisitionData, y, a_y); SIRF_DYNAMIC_CAST(const STIRAcquisitionData, b, a_b); - float a = *static_cast(ptr_a); int n = get_max_segment_num(); int nx = x.get_max_segment_num(); int ny = y.get_max_segment_num(); @@ -446,8 +446,8 @@ STIRImageData::write(const std::string &filename, const std::string &format_file format_sptr->write_to_file(filename, image); } -void -STIRImageData::sum(void* ptr) const +float +STIRImageData::sum() const { #if defined(_MSC_VER) && _MSC_VER < 1900 Image3DF::const_full_iterator iter; @@ -458,12 +458,13 @@ STIRImageData::sum(void* ptr) const double s = 0.0; for (iter = data().begin_all(); iter != data().end_all(); iter++) s += *iter; - float* ptr_s = static_cast(ptr); - *ptr_s = (float)s; + return (float)s; + //float* ptr_s = static_cast(ptr); + //*ptr_s = (float)s; } -void -STIRImageData::max(void* ptr) const +float +STIRImageData::max() const { #if defined(_MSC_VER) && _MSC_VER < 1900 Image3DF::const_full_iterator iter; @@ -474,12 +475,13 @@ STIRImageData::max(void* ptr) const float s = 0.0; for (iter = data().begin_all(); iter != data().end_all(); iter++) s = std::max(s, *iter); - float* ptr_s = static_cast(ptr); - *ptr_s = (float)s; + return s; + //float* ptr_s = static_cast(ptr); + //*ptr_s = (float)s; } -void -STIRImageData::dot(const DataContainer& a_x, void* ptr) const +float +STIRImageData::dot(const DataContainer& a_x) const { SIRF_DYNAMIC_CAST(const STIRImageData, x, a_x); #if defined(_MSC_VER) && _MSC_VER < 1900 @@ -497,26 +499,25 @@ STIRImageData::dot(const DataContainer& a_x, void* ptr) const double t = *iter; s += t * (*iter_x); } - float* ptr_s = static_cast(ptr); - *ptr_s = (float)s; + return (float)s; + //float* ptr_s = static_cast(ptr); + //*ptr_s = (float)s; } void STIRImageData::axpby( -const void* ptr_a, const DataContainer& a_x, -const void* ptr_b, const DataContainer& a_y) + float a, const DataContainer& a_x, + float b, const DataContainer& a_y) { //add deprecation warning - STIRImageData::xapyb(a_x, ptr_a, a_y, ptr_b); + STIRImageData::xapyb(a_x, a, a_y, b); } void STIRImageData::xapyb( -const DataContainer& a_x, const void* ptr_a, -const DataContainer& a_y, const void* ptr_b) + const DataContainer& a_x, float a, + const DataContainer& a_y, float b) { - float a = *static_cast(ptr_a); - float b = *static_cast(ptr_b); SIRF_DYNAMIC_CAST(const STIRImageData, x, a_x); SIRF_DYNAMIC_CAST(const STIRImageData, y, a_y); #if defined(_MSC_VER) && _MSC_VER < 1900 @@ -536,16 +537,15 @@ const DataContainer& a_y, const void* ptr_b) iter_x = x.data().begin_all(), iter_y = y.data().begin_all(); iter != data().end_all() && iter_x != x.data().end_all() && iter_y != y.data().end_all(); - iter++, iter_x++, iter_y++) + iter++, iter_x++, iter_y++) *iter = a * (*iter_x) + b * (*iter_y); } void STIRImageData::xapyb( -const DataContainer& a_x, const void* ptr_a, -const DataContainer& a_y, const DataContainer& a_b) + const DataContainer& a_x, float a, + const DataContainer& a_y, const DataContainer& a_b) { - float a = *static_cast(ptr_a); SIRF_DYNAMIC_CAST(const STIRImageData, b, a_b); SIRF_DYNAMIC_CAST(const STIRImageData, x, a_x); SIRF_DYNAMIC_CAST(const STIRImageData, y, a_y); diff --git a/src/xSTIR/cSTIR/stir_x.cpp b/src/xSTIR/cSTIR/stir_x.cpp index 9e4b18aac..80bbedd31 100644 --- a/src/xSTIR/cSTIR/stir_x.cpp +++ b/src/xSTIR/cSTIR/stir_x.cpp @@ -515,8 +515,8 @@ PETAcquisitionModel::forward(STIRAcquisitionData& ad, const STIRImageData& image if (sptr_add_.get() && !do_linear_only) { if (stir::Verbosity::get() > 1) std::cout << "additive term added..."; - ad.axpby(&one, ad, &one, *sptr_add_); - //ad.axpby(1.0, ad, 1.0, *sptr_add_); + //ad.axpby(&one, ad, &one, *sptr_add_); + ad.axpby(1.0, ad, 1.0, *sptr_add_); if (stir::Verbosity::get() > 1) std::cout << "ok\n"; } else @@ -533,8 +533,8 @@ PETAcquisitionModel::forward(STIRAcquisitionData& ad, const STIRImageData& image if (sptr_background_.get() && !do_linear_only) { if (stir::Verbosity::get() > 1) std::cout << "background term added..."; - ad.axpby(&one, ad, &one, *sptr_background_); - //ad.axpby(1.0, ad, 1.0, *sptr_background_); + //ad.axpby(&one, ad, &one, *sptr_background_); + ad.axpby(1.0, ad, 1.0, *sptr_background_); if (stir::Verbosity::get() > 1) std::cout << "ok\n"; } else diff --git a/src/xSTIR/cSTIR/tests/test1.cpp b/src/xSTIR/cSTIR/tests/test1.cpp index d55093979..29841ab65 100644 --- a/src/xSTIR/cSTIR/tests/test1.cpp +++ b/src/xSTIR/cSTIR/tests/test1.cpp @@ -245,7 +245,7 @@ int test1() float alpha = 1.0 / sim_norm; float beta = -alpha; acq_diff.axpby - (&alpha, sim_data, &beta, acq_data); + (alpha, sim_data, beta, acq_data); std::cout << "relative data difference: " << acq_diff.norm() << std::endl; // backproject the simulated data @@ -264,7 +264,7 @@ int test1() alpha = 1.0 / im_norm; beta = -1.0 / bd_norm; img_diff.axpby - (&alpha, image_data, &beta, back_data); + (alpha, image_data, beta, back_data); std::cout << "relative images difference: " << img_diff.norm() << std::endl; // compute the norm of the linear part of the acquisition model