From 4dfdcba1cbe9686d69458d6ae624e64c43fd1d30 Mon Sep 17 00:00:00 2001 From: Evgueni Ovtchinnikov Date: Fri, 26 May 2023 14:14:11 +0100 Subject: [PATCH] implemented safer C++ data type handling, fixes #1193 --- .../cReg/NiftiBasedRegistration.cpp | 2 +- src/Registration/cReg/NiftiImageData.cpp | 71 +- .../cReg/NiftiImageData3DTensor.cpp | 4 +- src/Registration/cReg/NiftyResampler.cpp | 24 +- src/Registration/cReg/Registration.cpp | 6 +- src/Registration/cReg/Resample.cpp | 8 +- src/Registration/cReg/cReg.cpp | 40 +- src/Registration/cReg/cReg_p.cpp | 18 +- .../include/sirf/Reg/NiftiBasedRegistration.h | 2 +- .../cReg/include/sirf/Reg/NiftiImageData.h | 1034 +++++++++-------- .../cReg/include/sirf/Reg/NiftyRegistration.h | 8 +- .../cReg/include/sirf/Reg/NiftyResampler.h | 8 +- .../cReg/include/sirf/Reg/Registration.h | 16 +- .../cReg/include/sirf/Reg/Resample.h | 27 +- .../cSyn/include/sirf/Syn/utilities.h | 10 +- src/Synergistic/cSyn/tests/test_conv_img.cpp | 2 +- src/Synergistic/cSyn/utilities.cpp | 12 +- src/Synergistic/sirf_convert_image_type.cpp | 18 +- src/Synergistic/sirf_registration.cpp | 10 +- src/Synergistic/sirf_resample.cpp | 16 +- src/Synergistic/tests/test_cSynergistic.cpp | 10 +- src/common/ImageData.cpp | 32 +- src/common/SIRF.py | 12 +- src/common/csirf.cpp | 350 +++--- src/common/include/sirf/Syn/utilities.h | 10 +- .../include/sirf/common/DataContainer.h | 244 +++- src/common/include/sirf/common/ImageData.h | 28 +- src/common/include/sirf/common/csirf.h | 14 +- .../cGadgetron/gadgetron_data_containers.cpp | 114 +- .../Gadgetron/gadgetron_data_containers.h | 58 +- src/xGadgetron/cGadgetron/tests/mrtests.cpp | 8 +- src/xSTIR/cSTIR/cstir.cpp | 2 +- .../include/sirf/STIR/stir_data_containers.h | 113 +- src/xSTIR/cSTIR/stir_data_containers.cpp | 54 +- 34 files changed, 1397 insertions(+), 988 deletions(-) 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..4301679b6 100644 --- a/src/Registration/cReg/NiftiImageData.cpp +++ b/src/Registration/cReg/NiftiImageData.cpp @@ -64,24 +64,24 @@ 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::operator=(const ImageData& to_copy) { // Check for self-assignment if (this != &to_copy) { @@ -213,10 +213,11 @@ 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) { + auto in_sptr = std::dynamic_pointer_cast >(dc_sptr); // Create image from input - out_sptr = std::make_shared >(*in_sptr); + out_sptr = std::make_shared >(*dc_sptr); auto &it_in = in_sptr->begin(); auto &it_out = out_sptr->begin(); @@ -225,13 +226,14 @@ 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"; // Create image from input - out_sptr = std::make_shared >(*in_sptr); + out_sptr = std::make_shared >(*dc_sptr); auto &it_in = in_sptr->begin(); auto &it_out = out_sptr->begin(); @@ -240,7 +242,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); @@ -1670,10 +1672,11 @@ 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); +// dataType s; +// this->dot(other, &s); +// return s; +//// return std::inner_product(this->begin(),this->end(),other.begin(),dataType(0)); } template @@ -1738,29 +1741,32 @@ 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; + //float* ptr_s = static_cast(ptr); + //*ptr_s = (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; + //float* ptr_s = static_cast(ptr); + //*ptr_s = (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,8 +1774,9 @@ void NiftiImageData::max(void* ptr) const if (si > s) s = si; } - float* ptr_s = static_cast(ptr); - *ptr_s = s; + return s; + //float* ptr_s = static_cast(ptr); + //*ptr_s = s; } template @@ -1919,56 +1926,56 @@ template void NiftiImageData::multiply (const DataContainer& a_x, const DataContainer& a_y) { - binary_op(a_x, a_y, DataContainer::product); + binary_op(a_x, a_y, DataContainerTempl::product); } template void NiftiImageData::multiply (const DataContainer& a_x, const void* a_y) { - semibinary_op(a_x, a_y, DataContainer::product); + semibinary_op(a_x, a_y, DataContainerTempl::product); } template void NiftiImageData::add (const DataContainer& a_x, const void* a_y) { - semibinary_op(a_x, a_y, DataContainer::sum); + semibinary_op(a_x, a_y, DataContainerTempl::sum); } template void NiftiImageData::divide (const DataContainer& a_x, const DataContainer& a_y) { - binary_op(a_x, a_y, DataContainer::ratio); + binary_op(a_x, a_y, DataContainerTempl::ratio); } template void NiftiImageData::maximum (const DataContainer& a_x, const DataContainer& a_y) { - binary_op(a_x, a_y, DataContainer::maximum); + binary_op(a_x, a_y, DataContainerTempl::maximum); } template void NiftiImageData::maximum (const DataContainer& a_x, const void* a_y) { - semibinary_op(a_x, a_y, DataContainer::maximum); + semibinary_op(a_x, a_y, DataContainerTempl::maximum); } template void NiftiImageData::minimum (const DataContainer& a_x, const DataContainer& a_y) { - binary_op(a_x, a_y, DataContainer::minimum); + binary_op(a_x, a_y, DataContainerTempl::minimum); } template void NiftiImageData::minimum (const DataContainer& a_x, const void* a_y) { - semibinary_op(a_x, a_y, DataContainer::minimum); + semibinary_op(a_x, a_y, DataContainerTempl::minimum); } template @@ -2006,13 +2013,13 @@ void NiftiImageData::sqrt(const DataContainer& a_x) template void NiftiImageData::sign(const DataContainer& a_x) { - unary_op(a_x, DataContainer::sign); + unary_op(a_x, DataContainerTempl::sign); } template void NiftiImageData::abs(const DataContainer& a_x) { - unary_op(a_x, DataContainer::abs); + unary_op(a_x, DataContainerTempl::abs); } template diff --git a/src/Registration/cReg/NiftiImageData3DTensor.cpp b/src/Registration/cReg/NiftiImageData3DTensor.cpp index 532cda4ed..9e3307add 100644 --- a/src/Registration/cReg/NiftiImageData3DTensor.cpp +++ b/src/Registration/cReg/NiftiImageData3DTensor.cpp @@ -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..48432ecdf 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,8 +229,9 @@ 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) { + std::shared_ptr > output_to_return_sptr = std::dynamic_pointer_cast >(out_to_return_sptr); // If output is only real, set that if (!output_to_return_sptr->is_complex()) output_to_return_sptr->fill(*resampled_niftis.real()); @@ -238,7 +240,7 @@ static void set_post_resample_outputs(std::shared_ptr &output_to_retu NumberType::Type output_num_type = (*output_to_return_sptr->begin()).get_typeID(); if (output_num_type != NumberType::CXFLOAT) throw std::runtime_error("NiftyResampler: Only complex type currently supported is complex float"); - ImageData::Iterator &it_out = output_to_return_sptr->begin(); + 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) { @@ -252,19 +254,19 @@ static void set_post_resample_outputs(std::shared_ptr &output_to_retu } 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 +297,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..1419832da 100644 --- a/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h +++ b/src/Registration/cReg/include/sirf/Reg/NiftiImageData.h @@ -41,621 +41,637 @@ 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 + 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 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 { + 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) + { + 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); - /// Filename constructor - NiftiImageData(const std::string &filename); + /// Assignment + 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); + //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); + //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); + //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; } + /// below all void* are actually float* + virtual float norm() const; + 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); + 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 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; + }; +} \ No newline at end of file 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..ddd2da487 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,19 +51,20 @@ 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); @@ -73,7 +75,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 +178,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..51d4356de 100644 --- a/src/Synergistic/tests/test_cSynergistic.cpp +++ b/src/Synergistic/tests/test_cSynergistic.cpp @@ -185,8 +185,8 @@ int main(int argc, char* argv[]) res_complex.set_floating_image(ismrmrd_im_sptr); 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::shared_ptr forward_cplx_sptr = res_complex.forward(ismrmrd_im_sptr); + std::shared_ptr adjoint_cplx_sptr = res_complex.adjoint(ismrmrd_im_sptr); // Get the output std::shared_ptr > forward_cplx_real_sptr, forward_cplx_imag_sptr, adjoint_cplx_real_sptr, adjoint_cplx_imag_sptr; @@ -297,10 +297,12 @@ 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::shared_ptr resampled_G2_sptr = res.forward(G2_sptr); std::cout << "\n reoriented back to original space:\n"; - resampled_G2_sptr->get_geom_info_sptr()->print_info(); + std::shared_ptr sptr_im = std::dynamic_pointer_cast(resampled_G2_sptr); + sptr_im->get_geom_info_sptr()->print_info(); + //resampled_G2_sptr->get_geom_info_sptr()->print_info(); if (NiftiImageData(*G1_sptr) != NiftiImageData(*resampled_G2_sptr)) throw std::runtime_error("GadgetronImagesVector::reorient test failed"); 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..fbc338b93 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,79 @@ 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); + auto& z = objectFromHandle >(h); + z.xapyb(x, ptr_a, y, ptr_b); +} + extern "C" void* cSIRF_axpby( @@ -185,11 +232,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 +242,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 +277,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 +286,32 @@ 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); + auto const& b = objectFromHandle >(ptr_b); + auto& z = objectFromHandle >(h); + z.xapyb(x, ptr_a, y, b); +} + extern "C" void* cSIRF_XapYB( @@ -258,11 +320,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 +329,35 @@ 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); + auto& z = objectFromHandle >(ptr_z); + z.add(x, ptr_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 +370,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 +412,34 @@ 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); + 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__); +} + extern "C" void* cSIRF_semibinary(const void* ptr_x, const void* ptr_y, const char* f) @@ -371,17 +447,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 +455,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 +481,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 +500,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 +543,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 +570,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** ptr; + SELECT_DATA_CASE(id.data_type(), geom_info_templ, ptr_im, ptr); + return *ptr; } 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..c02aa0c0a 100644 --- a/src/common/include/sirf/common/DataContainer.h +++ b/src/common/include/sirf/common/DataContainer.h @@ -37,11 +37,239 @@ which rely on the same features of the items. namespace sirf { - typedef std::map Dimensions; - class DataContainer { public: virtual ~DataContainer() {} + + 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 + { + // default value + return false; + } + + /// returns the size of data elements + virtual int bits() const + { + // default value + 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; + + 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()"); + } + }; + + 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 T sum() const = 0; + + /// calculates the value of this container's element with the largest real part + 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; + /// \c *this = the product \c x * y with scalar y + virtual void multiply + (const DataContainer& x, const void* ptr_y) = 0; + + /// \c *this = the sum \c x + y with scalar y + virtual void add + (const DataContainer& x, const void* ptr_y) = 0; + + /// \c *this = the elementwise ratio \c x / y + virtual void divide + (const DataContainer& x, const DataContainer& y) = 0; + + /// \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; + + /// \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; + + /// \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; + + /// \c *this = the elementwise \c exp(x) + virtual void exp(const DataContainer& x) = 0; + /// \c *this = the elementwise \c log(x) + virtual void log(const DataContainer& x) = 0; + /// \c *this = the elementwise \c sqrt(x) + virtual void sqrt(const DataContainer& x) = 0; + /// \c *this = the elementwise \c sign(x) + virtual void sign(const DataContainer& x) = 0; + /// \c *this = the elementwise \c abs(x) + virtual void abs(const DataContainer& x) = 0; + + /// \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 + 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( + const DataContainer& x, const DataContainer& a, + const DataContainer& y, const DataContainer& b) = 0; + + /// \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_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); + } + + static T product(T x, T y) + { + return x * y; + } + + static T ratio(T x, T y) + { + return x / y; + } + + static T inverse_ratio(T x, T y) + { + return y / x; + } + + static T sum(T x, T y) + { + return x + y; + } + + static T maximum(T x, T y) + { + return std::max(x, y); + } + static T maxabs(T x, T y) + { + return std::max(std::abs(x), std::abs(y)); + } + static T maxreal(T x, T y) + { + return std::real(x) > std::real(y) ? x : y; + } + + static T minimum(T x, T y) + { + return std::min(x, y); + } + static T minabs(T x, T y) + { + return std::min(std::abs(x), std::abs(y)); + } + 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) + { + return std::pow(x, y); + } + static std::complex exp(std::complex x) + { + return std::exp(x); + } + static std::complex log(std::complex x) + { + return std::log(x); + } + static std::complex sqrt(std::complex x) + { + return std::sqrt(x); + } + static T sign(T x) + { + return (std::real(x) > 0) - (std::real(x) < 0); + } + static T abs(T x) + { + return T(std::abs(x)); + } + }; + + typedef std::map Dimensions; + + // Old DataContainer - to be removed + class OldDataContainer { + public: + virtual ~OldDataContainer() {} //virtual DataContainer* new_data_container() const = 0; virtual ObjectHandle* new_data_container_handle() const = 0; virtual unsigned int items() const = 0; @@ -146,9 +374,9 @@ namespace sirf { return items() < 1; } - std::unique_ptr clone() const + std::unique_ptr clone() const { - return std::unique_ptr(this->clone_impl()); + return std::unique_ptr(this->clone_impl()); } /// overwrites this container's complex data with complex conjugate values @@ -158,11 +386,11 @@ namespace sirf { } /// returns unique pointer to the complex-conjugated copy of this container - std::unique_ptr conjugate() const + std::unique_ptr conjugate() const { - DataContainer* ptr = this->clone_impl(); + OldDataContainer* ptr = this->clone_impl(); ptr->conjugate(); - return std::unique_ptr(ptr); + return std::unique_ptr(ptr); } template @@ -248,7 +476,7 @@ namespace sirf { } protected: - virtual DataContainer* clone_impl() const = 0; + virtual OldDataContainer* clone_impl() const = 0; /// we assume data to be real, complex data containers must override this virtual void conjugate_impl() { 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..8e31ade92 100644 --- a/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp +++ b/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp @@ -417,105 +417,105 @@ 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 @@ -570,8 +570,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 +592,13 @@ MRAcquisitionData::dot(const DataContainer& dc, void* ptr) const i++; j++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = z; } -void -MRAcquisitionData::sum(void* ptr) const +complex_float_t +MRAcquisitionData::sum() const { int n = number(); complex_float_t z = 0; @@ -610,12 +611,13 @@ MRAcquisitionData::sum(void* ptr) const z += MRAcquisitionData::sum(a); i++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = z; } -void -MRAcquisitionData::max(void* ptr) const +complex_float_t +MRAcquisitionData::max() const { int n = number(); complex_float_t z = 0; @@ -632,8 +634,9 @@ MRAcquisitionData::max(void* ptr) const z = zi; i++; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = z; } void @@ -1416,8 +1419,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 +1429,13 @@ 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; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = 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 +1443,13 @@ 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; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = 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,8 +1460,9 @@ GadgetronImageData::max(void* ptr) const if (ri > r) z = zi; } - complex_float_t* ptr_z = static_cast(ptr); - *ptr_z = z; + return z; + //complex_float_t* ptr_z = static_cast(ptr); + //*ptr_z = z; } void @@ -1571,7 +1577,7 @@ 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); + binary_op(x, y, DataContainerTempl::product); } void @@ -1579,7 +1585,7 @@ 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); + semibinary_op(x, y, DataContainerTempl::product); } void @@ -1587,7 +1593,7 @@ 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); + semibinary_op(x, y, DataContainerTempl::sum); } void @@ -1595,7 +1601,7 @@ 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); + binary_op(x, y, DataContainerTempl::ratio); } void @@ -1605,7 +1611,7 @@ GadgetronImageData::maximum( { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::maxreal); + binary_op(x, y, DataContainerTempl::maxreal); } void @@ -1613,7 +1619,7 @@ 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); + semibinary_op(x, y, DataContainerTempl::maxreal); } void @@ -1623,7 +1629,7 @@ GadgetronImageData::minimum( { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::minreal); + binary_op(x, y, DataContainerTempl::minreal); } void @@ -1631,7 +1637,7 @@ 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); + semibinary_op(x, y, DataContainerTempl::minreal); } void @@ -1641,7 +1647,7 @@ GadgetronImageData::power( { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); SIRF_DYNAMIC_CAST(const GadgetronImageData, y, a_y); - binary_op(x, y, DataContainer::power); + binary_op(x, y, DataContainerTempl::power); } void @@ -1649,42 +1655,42 @@ 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); + semibinary_op(x, y, DataContainerTempl::power); } void GadgetronImageData::exp(const DataContainer& a_x) { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::exp); + unary_op(x, DataContainerTempl::exp); } void GadgetronImageData::log(const DataContainer& a_x) { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::log); + unary_op(x, DataContainerTempl::log); } void GadgetronImageData::sqrt(const DataContainer& a_x) { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::sqrt); + unary_op(x, DataContainerTempl::sqrt); } void GadgetronImageData::sign(const DataContainer& a_x) { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::sign); + unary_op(x, DataContainerTempl::sign); } void GadgetronImageData::abs(const DataContainer& a_x) { SIRF_DYNAMIC_CAST(const GadgetronImageData, x, a_x); - unary_op(x, DataContainer::abs); + unary_op(x, DataContainerTempl::abs); } float 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..1a5288dd5 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,13 @@ namespace sirf { \brief Abstract MR acquisition data container class. */ - class MRAcquisitionData : public DataContainer { + class MRAcquisitionData : public DataContainerTempl { public: + virtual std::string data_type() const + { + return std::string("complex_float"); + } + // static methods // ISMRMRD acquisitions algebra: acquisitions viewed as vectors of @@ -538,15 +543,17 @@ 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 complex_float_t sum() const; + virtual complex_float_t max() const; + virtual complex_float_t dot(const DataContainer& dc) const; + //virtual void sum(void* ptr) const; + //virtual void max(void* ptr) const; + //complex_float_t dot(const DataContainer& a_x) + //{ + // complex_float_t z; + // dot(a_x, &z); + // return z; + //} virtual void axpby( const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y); @@ -716,7 +723,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,7 +789,7 @@ namespace sirf { */ - class ISMRMRDImageData : public ImageData { + class ISMRMRDImageData : public ImageData { public: //ISMRMRDImageData(ISMRMRDImageData& id, const char* attr, //const char* target); //does not build, have to be in the derived class @@ -886,9 +893,12 @@ 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 sum(void* ptr) const; + //virtual void max(void* ptr) const; + //virtual void dot(const DataContainer& dc, void* ptr) 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); @@ -950,12 +960,12 @@ namespace sirf { 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; - } + //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) @@ -1083,6 +1093,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 +1276,7 @@ namespace sirf { } virtual unsigned int items() const { - return (unsigned int)images_.size(); + return (unsigned int)this->images_.size(); } virtual unsigned int number() const { 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..b16cb88f7 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,29 +250,31 @@ 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 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& a_x) const; + //{ + // float s; + // dot(a_x, &s); + // return s; + //} virtual void axpby( const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y); @@ -299,32 +305,32 @@ namespace sirf { } virtual void sign(const DataContainer& x) { - unary_op(x, DataContainer::sign); + unary_op(x, DataContainerTempl::sign); } virtual void multiply(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::product); + semibinary_op(x, y, DataContainerTempl::product); } virtual void add(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::sum); + semibinary_op(x, y, DataContainerTempl::sum); } virtual void divide(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::ratio); + semibinary_op(x, y, DataContainerTempl::ratio); } virtual void maximum(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::maximum); + semibinary_op(x, y, DataContainerTempl::maximum); } virtual void minimum(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::minimum); + semibinary_op(x, y, DataContainerTempl::minimum); } virtual void power(const DataContainer& x, const void* ptr_y) { @@ -333,19 +339,19 @@ namespace sirf { } virtual void multiply(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::product); + binary_op(x, y, DataContainerTempl::product); } virtual void divide(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::ratio); + binary_op(x, y, DataContainerTempl::ratio); } virtual void maximum(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::maximum); + binary_op(x, y, DataContainerTempl::maximum); } virtual void minimum(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::minimum); + binary_op(x, y, DataContainerTempl::minimum); } virtual void power(const DataContainer& x, const DataContainer& y) { @@ -761,7 +767,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 +775,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 +783,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) { @@ -852,8 +858,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 { @@ -1049,9 +1059,12 @@ namespace sirf { 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 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); @@ -1082,32 +1095,32 @@ namespace sirf { } virtual void sign(const DataContainer& x) { - unary_op(x, DataContainer::sign); + unary_op(x, DataContainerTempl::sign); } virtual void multiply(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::product); + semibinary_op(x, y, DataContainerTempl::product); } virtual void add(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::sum); + semibinary_op(x, y, DataContainerTempl::sum); } virtual void divide(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::ratio); + semibinary_op(x, y, DataContainerTempl::ratio); } virtual void maximum(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::maximum); + semibinary_op(x, y, DataContainerTempl::maximum); } virtual void minimum(const DataContainer& x, const void* ptr_y) { float y = *static_cast(ptr_y); - semibinary_op(x, y, DataContainer::minimum); + semibinary_op(x, y, DataContainerTempl::minimum); } virtual void power(const DataContainer& x, const void* ptr_y) { @@ -1116,19 +1129,19 @@ namespace sirf { } virtual void multiply(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::product); + binary_op(x, y, DataContainerTempl::product); } virtual void divide(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::ratio); + binary_op(x, y, DataContainerTempl::ratio); } virtual void maximum(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::maximum); + binary_op(x, y, DataContainerTempl::maximum); } virtual void minimum(const DataContainer& x, const DataContainer& y) { - binary_op(x, y, DataContainer::minimum); + binary_op(x, y, DataContainerTempl::minimum); } virtual void power(const DataContainer& x, const DataContainer& y) { @@ -1169,12 +1182,12 @@ namespace sirf { _data->fill(v); } void scale(float s); - float dot(const DataContainer& a_x) const - { - float s; - dot(a_x, &s); - return 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) diff --git a/src/xSTIR/cSTIR/stir_data_containers.cpp b/src/xSTIR/cSTIR/stir_data_containers.cpp index 5211a10c6..c7127b13c 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,8 +126,9 @@ 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 @@ -446,8 +449,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 +461,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 +478,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,8 +502,9 @@ 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