[ismrmrd] 89/177: working on ndarray c++ api and utilities.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:02:06 UTC 2015
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to annotated tag v1.1.0.beta.1
in repository ismrmrd.
commit 87ee26f97a4cf343236d817c58ae3882ab1516ab
Author: Souheil Inati <souheil.inati at nih.gov>
Date: Thu Sep 25 09:01:48 2014 -0400
working on ndarray c++ api and utilities.
---
include/ismrmrd/dataset.h | 7 ++-
include/ismrmrd/ismrmrd.h | 23 ++++---
libsrc/dataset.cpp | 35 +++++++++--
libsrc/ismrmrd.cpp | 94 +++++++++++++++++++++++-----
utilities/generate_cartesian_shepp_logan.cpp | 67 +++++++++-----------
utilities/ismrmrd_fftw.h | 29 +++++----
utilities/ismrmrd_ndarray.h | 91 ---------------------------
utilities/ismrmrd_phantom.cpp | 34 +++++-----
utilities/ismrmrd_phantom.h | 16 ++---
9 files changed, 191 insertions(+), 205 deletions(-)
diff --git a/include/ismrmrd/dataset.h b/include/ismrmrd/dataset.h
index 2f238a2..38faf82 100644
--- a/include/ismrmrd/dataset.h
+++ b/include/ismrmrd/dataset.h
@@ -172,18 +172,19 @@ public:
uint32_t getNumberOfAcquisitions();
// Images
int appendImage(const std::string &var, const ISMRMRD_BlockModes blockmode, const Image &im);
+ int appendImage(const std::string &var, const ISMRMRD_BlockModes blockmode, const ISMRMRD_Image *im);
int readImage(const std::string &var, uint32_t index, Image &im);
uint32_t getNumberOfImages(const std::string &var);
// NDArrays
- int appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray &arr);
- int readNDArray(const std::string &var, uint32_t index, NDArray &arr);
+ template <typename T> int appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<T> &arr);
+ int appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const ISMRMRD_NDArray *arr);
+ template <typename T> int readNDArray(const std::string &var, uint32_t index, NDArray<T> &arr);
uint32_t getNumberOfNDArrays(const std::string &var);
protected:
ISMRMRD_Dataset dset_;
};
-
} /* ISMRMRD namespace */
#endif
diff --git a/include/ismrmrd/ismrmrd.h b/include/ismrmrd/ismrmrd.h
index dccfacf..ac4624e 100644
--- a/include/ismrmrd/ismrmrd.h
+++ b/include/ismrmrd/ismrmrd.h
@@ -339,8 +339,8 @@ typedef struct ISMRMRD_NDArray {
uint16_t version; /**< First unsigned int indicates the version */
uint16_t data_type; /**< e.g. unsigned short, float, complex float, etc. */
uint16_t ndim; /**< Number of dimensions */
- uint16_t dims[ISMRMRD_NDARRAY_MAXDIM]; /**< Dimensions */
- void *data; /**< Pointer to data */
+ size_t dims[ISMRMRD_NDARRAY_MAXDIM]; /**< Dimensions */
+ void *data; /**< Pointer to data */
} ISMRMRD_NDArray;
/** @addtogroup capi
@@ -586,24 +586,29 @@ public:
};
/// N-Dimensional array type
-class EXPORTISMRMRD NDArray: protected ISMRMRD_NDArray {
+template <typename T> class EXPORTISMRMRD NDArray: protected ISMRMRD_NDArray {
+ friend class Dataset;
public:
// Constructors, destructor and copy
NDArray();
- NDArray(const ISMRMRD_DataTypes dtype, const std::vector<uint16_t> dimvec);
- NDArray(const NDArray &other);
+ NDArray(const std::vector<size_t> dimvec);
+ NDArray(const NDArray<T> &other);
~NDArray();
- NDArray & operator= (const NDArray &other);
+ NDArray<T> & operator= (const NDArray<T> &other);
// Accessors and mutators
const uint16_t getVersion();
const ISMRMRD_DataTypes getDataType();
const uint16_t getNDim();
- const uint16_t (&getDims())[ISMRMRD_NDARRAY_MAXDIM];
- int setProperties(const ISMRMRD_DataTypes dtype, const std::vector<uint16_t> dimvec);
- void * getData();
+ const size_t (&getDims())[ISMRMRD_NDARRAY_MAXDIM];
+ const size_t getDataSize();
+ int resize(const std::vector<size_t> dimvec);
+ const size_t getNumberOfElements();
+ T * getData();
+
};
+
/** @} */
} // namespace ISMRMRD
diff --git a/libsrc/dataset.cpp b/libsrc/dataset.cpp
index 85320b8..ae6d4a7 100644
--- a/libsrc/dataset.cpp
+++ b/libsrc/dataset.cpp
@@ -77,6 +77,15 @@ int Dataset::appendImage(const std::string &var, const ISMRMRD_BlockModes blockm
return status;
}
+int Dataset::appendImage(const std::string &var, const ISMRMRD_BlockModes blockmode, const ISMRMRD_Image *im)
+{
+ int status = ismrmrd_append_image(&dset_, var.c_str(), blockmode, im);
+ if (status != ISMRMRD_NOERROR) {
+ //TODO throw an exception
+ }
+ return status;
+}
+
int Dataset::readImage(const std::string &var, uint32_t index, Image &im) {
int status = ismrmrd_read_image(&dset_, var.c_str(), index, reinterpret_cast<ISMRMRD_Image*>(&im));
if (status != ISMRMRD_NOERROR) {
@@ -93,17 +102,35 @@ uint32_t Dataset::getNumberOfImages(const std::string &var)
// NDArrays
-int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray &arr)
+template <typename T> int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<T> &arr)
+{
+ int status = ismrmrd_append_array(&dset_, var.c_str(), blockmode, static_cast<const ISMRMRD_NDArray*>(&arr));
+ if (status != ISMRMRD_NOERROR) {
+ //TODO throw an exception
+ }
+ return status;
+}
+// Specific instantiations
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<uint16_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<int16_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<uint32_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<int32_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<float_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<double_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<complex_float_t> &arr);
+template int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const NDArray<complex_double_t> &arr);
+
+int Dataset::appendNDArray(const std::string &var, const ISMRMRD_BlockModes blockmode, const ISMRMRD_NDArray *arr)
{
- int status = ismrmrd_append_array(&dset_, var.c_str(), blockmode, reinterpret_cast<const ISMRMRD_NDArray*>(&arr));
+ int status = ismrmrd_append_array(&dset_, var.c_str(), blockmode, arr);
if (status != ISMRMRD_NOERROR) {
//TODO throw an exception
}
return status;
}
-int Dataset::readNDArray(const std::string &var, uint32_t index, NDArray &arr) {
- int status = ismrmrd_read_array(&dset_, var.c_str(), index, reinterpret_cast<ISMRMRD_NDArray*>(&arr));
+template <typename T> int Dataset::readNDArray(const std::string &var, uint32_t index, NDArray<T> &arr) {
+ int status = ismrmrd_read_array(&dset_, var.c_str(), index, static_cast<ISMRMRD_NDArray*>(&arr));
if (status != ISMRMRD_NOERROR) {
//TODO throw an exception
}
diff --git a/libsrc/ismrmrd.cpp b/libsrc/ismrmrd.cpp
index 29485e3..6b0a6b9 100644
--- a/libsrc/ismrmrd.cpp
+++ b/libsrc/ismrmrd.cpp
@@ -5,6 +5,48 @@
namespace ISMRMRD {
+// Internal function to control the allowed data types for NDArrays
+template <typename T> ISMRMRD_DataTypes get_data_type();
+template <> ISMRMRD_DataTypes get_data_type<uint16_t>()
+{
+ return ISMRMRD_USHORT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<int16_t>()
+{
+ return ISMRMRD_SHORT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<uint32_t>()
+{
+ return ISMRMRD_UINT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<int32_t>()
+{
+ return ISMRMRD_INT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<float>()
+{
+ return ISMRMRD_FLOAT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<double>()
+{
+ return ISMRMRD_DOUBLE;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<complex_float_t>()
+{
+ return ISMRMRD_CXFLOAT;
+}
+
+template <> inline ISMRMRD_DataTypes get_data_type<complex_double_t>()
+{
+ return ISMRMRD_CXDOUBLE;
+}
+
//
// AcquisitionHeader class implementation
//
@@ -447,29 +489,31 @@ void Image::clearAllFlags() {
//
// Array class Implementation
//
-NDArray::NDArray()
+template <typename T> NDArray<T>::NDArray()
{
ismrmrd_init_ndarray(this);
+ data_type = get_data_type<T>();
}
-NDArray::NDArray(const ISMRMRD_DataTypes dtype, const std::vector<uint16_t> dimvec)
+template <typename T> NDArray<T>::NDArray(const std::vector<size_t> dimvec)
{
ismrmrd_init_ndarray(this);
- setProperties(dtype, dimvec);
+ data_type = get_data_type<T>();
+ resize(dimvec);
}
-NDArray::NDArray(const NDArray &other)
+template <typename T> NDArray<T>::NDArray(const NDArray<T> &other)
{
ismrmrd_init_ndarray(this);
ismrmrd_copy_ndarray(this, &other);
}
-NDArray::~NDArray()
+template <typename T> NDArray<T>::~NDArray()
{
ismrmrd_cleanup_ndarray(this);
}
-NDArray & NDArray::operator= (const NDArray &other)
+template <typename T> NDArray<T> & NDArray<T>::operator= (const NDArray<T> &other)
{
// Assignment makes a copy
if (this != &other )
@@ -480,28 +524,27 @@ NDArray & NDArray::operator= (const NDArray &other)
return *this;
}
-const uint16_t NDArray::getVersion() {
+template <typename T> const uint16_t NDArray<T>::getVersion() {
return version;
};
-const ISMRMRD_DataTypes NDArray::getDataType() {
+template <typename T> const ISMRMRD_DataTypes NDArray<T>::getDataType() {
return static_cast<ISMRMRD_DataTypes>( data_type );
}
-const uint16_t NDArray::getNDim() {
+template <typename T> const uint16_t NDArray<T>::getNDim() {
return ndim;
};
-const uint16_t (&NDArray::getDims())[ISMRMRD_NDARRAY_MAXDIM] {
+template <typename T> const size_t (&NDArray<T>::getDims())[ISMRMRD_NDARRAY_MAXDIM] {
return dims;
};
-int NDArray::setProperties(const ISMRMRD_DataTypes dtype, const std::vector<uint16_t> dimvec) {
+template <typename T> int NDArray<T>::resize(const std::vector<size_t> dimvec) {
if (dimvec.size() > ISMRMRD_NDARRAY_MAXDIM) {
// TODO throw exception
return ISMRMRD_MEMORYERROR;
}
- data_type = dtype;
ndim = dimvec.size();
for (int n=0; n<ndim; n++) {
dims[n] = dimvec[n];
@@ -510,8 +553,31 @@ int NDArray::setProperties(const ISMRMRD_DataTypes dtype, const std::vector<uint
return status;
}
-void * NDArray::getData() {
- return data;
+template <typename T> T * NDArray<T>::getData() {
+ return static_cast<T*>(data);
+}
+
+template <typename T> const size_t NDArray<T>::getDataSize() {
+ return ismrmrd_size_of_ndarray_data(this);
+}
+
+template <typename T> const size_t NDArray<T>::getNumberOfElements() {
+ size_t num = 1;
+ for (int n = 0; n < ndim; n++) {
+ num *= dims[n];
+ }
+ return num;
}
+
+// Specializations
+template class NDArray<uint16_t>;
+template class NDArray<int16_t>;
+template class NDArray<uint32_t>;
+template class NDArray<int32_t>;
+template class NDArray<float_t>;
+template class NDArray<double_t>;
+template class NDArray<complex_float_t>;
+template class NDArray<complex_double_t>;
+
} // namespace ISMRMRD
diff --git a/utilities/generate_cartesian_shepp_logan.cpp b/utilities/generate_cartesian_shepp_logan.cpp
index a290ef1..975bdb5 100644
--- a/utilities/generate_cartesian_shepp_logan.cpp
+++ b/utilities/generate_cartesian_shepp_logan.cpp
@@ -66,16 +66,16 @@ int main(int argc, char** argv)
std::cout << "Generating Cartesian Shepp Logan Phantom!!!" << std::endl;
std::cout << "Acceleration: " << acc_factor << std::endl;
- boost::shared_ptr<NDArrayContainer<std::complex<float> > > phantom = shepp_logan_phantom(matrix_size);
- boost::shared_ptr<NDArrayContainer<std::complex<float> > > coils = generate_birdcage_sensititivies(matrix_size, ncoils, 1.5);
+ boost::shared_ptr<NDArray<complex_float_t> > phantom = shepp_logan_phantom(matrix_size);
+ boost::shared_ptr<NDArray<complex_float_t> > coils = generate_birdcage_sensititivies(matrix_size, ncoils, 1.5);
- std::vector<unsigned int> dims;
+ std::vector<size_t> dims;
dims.push_back(matrix_size*ros); //oversampling in the readout direction
dims.push_back(matrix_size);
dims.push_back(ncoils);
- NDArrayContainer<std::complex<float> > coil_images(dims);
- coil_images.data_ = std::complex<float>(0.0,0.0);
+ NDArray<complex_float_t> coil_images(dims);
+ memset(coil_images.getData(), 0, coil_images.getDataSize());
for (unsigned int c = 0; c < ncoils; c++) {
for (unsigned int y = 0; y < matrix_size; y++) {
@@ -83,13 +83,12 @@ int main(int argc, char** argv)
size_t out_index = c*matrix_size*matrix_size*ros + y*matrix_size*ros + ((matrix_size*ros-matrix_size)>>1) + x;
size_t cindex = c*matrix_size*matrix_size + y*matrix_size + x;
size_t iindex = y*matrix_size + x;
- coil_images.data_[out_index] = phantom->data_[iindex] * coils->data_[cindex];
+ coil_images.getData()[out_index] = phantom->getData()[iindex] * coils->getData()[cindex];
}
}
}
-
- //Let's append the data to the file
+ //Let's append the data to the file
//Create if needed
Dataset d(outfile.c_str(),dataset.c_str(), true);
Acquisition acq;
@@ -110,10 +109,18 @@ int main(int argc, char** argv)
d.appendAcquisition(acq);
}
- for (unsigned int r = 0; r < repetitions; r++) {
+ //NDArray<complex_float_t> cm(dims);
+ std::cout << "FFT number 0" << std::endl;
+ for (unsigned int r = 0; r < repetitions; r++) {
for (unsigned int a = 0; a < acc_factor; a++) {
- NDArrayContainer<std::complex<float> > cm = coil_images;
+ //for (size_t n=0; n<coil_images.getNumberOfElements(); n++) {
+ // cm.getData()[n] = coil_images.getData()[n];
+ //}
+ NDArray<complex_float_t> cm = coil_images;
fft2c(cm);
+ std::cout << "FFT number 1" << std::endl;
+
+
add_noise(cm,noise_level);
for (int i = a; i < matrix_size; i+=acc_factor) {
acq.clearAllFlags();
@@ -131,7 +138,7 @@ int main(int argc, char** argv)
acq.sample_time_us() = 5.0;
for (unsigned int c = 0; c < ncoils; c++) {
memcpy(&(acq.getData()[c*readout]),
- &(cm.data_[c*matrix_size*readout + i*readout]),
+ &(cm.getData()[c*matrix_size*readout + i*readout]),
sizeof(complex_float_t)*readout);
}
@@ -149,6 +156,9 @@ int main(int argc, char** argv)
}
}
+ std::cout << "Wrote acquisitions" << std::endl;
+
+
//Let's create a header, we will use the C++ classes in ismrmrd/xml.h
IsmrmrdHeader h;
h.experimentalConditions.H1resonanceFrequency_Hz = 63500000; //~1.5T
@@ -199,34 +209,13 @@ int main(int argc, char** argv)
//Write the header to the data file.
d.writeHeader(xml_header);
- // Convert the NDArrayContainer to an ISMRMRD::NDArray before appending
- // This requires a memcopy and is annoying
- // TODO fix this
- NDArray arr;
- std::vector<uint16_t> arrdims;
- arrdims.resize(phantom->ndims());
- for (int n=0; n<phantom->ndims(); n++) {
- arrdims[n] = phantom->dimensions_[n];
- }
- arr.setProperties(ISMRMRD_CXFLOAT, arrdims);
- memcpy(arr.getData(), &(phantom->data_[0]), phantom->elements()*sizeof(complex_float_t));
- d.appendNDArray("phantom", ISMRMRD_BLOCKMODE_ARRAY, arr);
-
- arrdims.resize(coils->ndims());
- for (int n=0; n<coils->ndims(); n++) {
- arrdims[n] = coils->dimensions_[n];
- }
- arr.setProperties(ISMRMRD_CXFLOAT, arrdims);
- memcpy(arr.getData(), &(coils->data_[0]), coils->elements()*sizeof(complex_float_t));
- d.appendNDArray("csm", ISMRMRD_BLOCKMODE_ARRAY, arr);
-
- arrdims.resize(coil_images.ndims());
- for (int n=0; n<coil_images.ndims(); n++) {
- arrdims[n] = coil_images.dimensions_[n];
- }
- arr.setProperties(ISMRMRD_CXFLOAT, arrdims);
- memcpy(arr.getData(), &(coil_images.data_[0]), coil_images.elements()*sizeof(complex_float_t));
- d.appendNDArray("coil_images", ISMRMRD_BLOCKMODE_ARRAY, arr);
+ std::cout << "Generating Cartesian Shepp Logan Phantom!!!" << std::endl;
+
+ //Write out some arrays for convenience
+ d.appendNDArray("phantom", ISMRMRD_BLOCKMODE_ARRAY, *phantom);
+ d.appendNDArray("csm", ISMRMRD_BLOCKMODE_ARRAY, *coils);
+ d.appendNDArray("coil_images", ISMRMRD_BLOCKMODE_ARRAY, coil_images);
+
return 0;
}
diff --git a/utilities/ismrmrd_fftw.h b/utilities/ismrmrd_fftw.h
index 94e68b1..45ac33c 100644
--- a/utilities/ismrmrd_fftw.h
+++ b/utilities/ismrmrd_fftw.h
@@ -6,7 +6,6 @@
*/
#include "fftw3.h"
-#include "ismrmrd_ndarray.h"
namespace ISMRMRD {
@@ -24,56 +23,56 @@ template<typename TI, typename TO> void circshift(TO *out, const TI *in, int xdi
#define fftshift(out, in, x, y) circshift(out, in, x, y, (x/2), (y/2))
-int fft2c(NDArrayContainer<std::complex<float> >& a, bool forward)
+int fft2c(NDArray<complex_float_t> &a, bool forward)
{
- if (a.dimensions_.size() < 2) {
+ if (a.getNDim() < 2) {
std::cout << "fft2c Error: input array must have at least two dimensions" << std::endl;
return -1;
}
- size_t elements = a.dimensions_[0]*a.dimensions_[1];
- size_t ffts = a.data_.size()/elements;
+ size_t elements = a.getDims()[0]*a.getDims()[1];
+ size_t ffts = a.getNumberOfElements()/elements;
//Array for transformation
- fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*a.data_.size());
+ fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*a.getNumberOfElements());
if (!tmp) {
std::cout << "Error allocating temporary storage for FFTW" << std::endl;
return -1;
}
-
for (size_t f = 0; f < ffts; f++) {
size_t index = f*elements;
- fftshift(reinterpret_cast<std::complex<float>*>(tmp),&a.data_[index],a.dimensions_[0],a.dimensions_[1]);
+ fftshift(reinterpret_cast<std::complex<float>*>(tmp),&a.getData()[index],a.getDims()[0],a.getDims()[1]);
//Create the FFTW plan
fftwf_plan p;
if (forward) {
- p = fftwf_plan_dft_2d(a.dimensions_[1], a.dimensions_[0], tmp,tmp, FFTW_FORWARD, FFTW_ESTIMATE);
+ p = fftwf_plan_dft_2d(a.getDims()[1], a.getDims()[0], tmp,tmp, FFTW_FORWARD, FFTW_ESTIMATE);
} else {
- p = fftwf_plan_dft_2d(a.dimensions_[1], a.dimensions_[0], tmp,tmp, FFTW_BACKWARD, FFTW_ESTIMATE);
+ p = fftwf_plan_dft_2d(a.getDims()[1], a.getDims()[0], tmp,tmp, FFTW_BACKWARD, FFTW_ESTIMATE);
}
fftwf_execute(p);
- fftshift(&a.data_[index],reinterpret_cast<std::complex<float>*>(tmp),a.dimensions_[0],a.dimensions_[1]);
+ fftshift(&a.getData()[index],reinterpret_cast<std::complex<float>*>(tmp),a.getDims()[0],a.getDims()[1]);
//Clean up.
fftwf_destroy_plan(p);
}
std::complex<float> scale(std::sqrt(1.0f*elements),0.0);
- a.data_ /= scale;
-
+ for (size_t n=0; n<a.getNumberOfElements(); n++) {
+ a.getData()[n] /= scale;
+ }
fftwf_free(tmp);
return 0;
}
-int fft2c(NDArrayContainer<std::complex<float> >& a) {
+int fft2c(NDArray<complex_float_t> &a) {
return fft2c(a,true);
}
-int ifft2c(NDArrayContainer<std::complex<float> >& a) {
+int ifft2c(NDArray<complex_float_t> &a) {
return fft2c(a,false);
}
diff --git a/utilities/ismrmrd_ndarray.h b/utilities/ismrmrd_ndarray.h
deleted file mode 100644
index f37c57f..0000000
--- a/utilities/ismrmrd_ndarray.h
+++ /dev/null
@@ -1,91 +0,0 @@
-/**
- * Container for generic multi-dimensional array.
- */
-
-#pragma once
-#ifndef ISMRMRD_NDARRAY_H
-#define ISMRMRD_NDARRAY_H
-
-#include <valarray>
-
-namespace ISMRMRD {
-
-template <typename T> class NDArrayContainer {
-
-public:
- NDArrayContainer() {}
-
- /**
- * @brief Construct with dimensions and data
- */
- NDArrayContainer(const std::vector<unsigned int>& dimensions) {
- dimensions_ = dimensions;
- data_.resize(elements());
- }
-
- /**
- * @brief Construct with dimensions and data
- */
- NDArrayContainer(const std::vector<unsigned int>& dimensions, const T* d) {
- dimensions_ = dimensions;
- data_.resize(elements());
- memcpy(&data_[0],d,sizeof(T)*elements());
- }
-
- /**
- * @brief Construct with dimensions and preset value
- */
- NDArrayContainer(const std::vector<unsigned int>& dimensions, const std::valarray<T>& t) {
- dimensions_ = dimensions;
- data_.resize(elements());
- data_ = t;
- }
-
- virtual ~NDArrayContainer() {}
-
- std::vector<unsigned int> dimensions_; /**< Array with dimensions of the array. First dimension is fastest moving in the array */
- std::valarray<T> data_; /**< The data itself. A vector is used here for easy memory management */
-
- size_t elements() const {
- if (dimensions_.size() == 0) {
- return 0;
- }
- size_t elements = 1;
- std::vector<unsigned int>::const_iterator it = dimensions_.begin();
- while (it != dimensions_.end()) {
- elements *= *(it++);
- }
- return elements;
- }
-
- T& operator[] (const size_t& p) {
- return data_[p];
- }
-
- T operator[] (const size_t& p) const {
- return data_[p];
- }
-
- void resize (const size_t& s) {
- data_.resize(s);
- }
-
- void resize (const size_t& s, const T& t) {
- data_.resize(s,t);
- }
-
- bool is_consistent() const {
- return (elements() == data_.size());
- }
-
- size_t ndims() const {
- size_t i = 1, n_dimensions = 1;
- for (; i < dimensions_.size(); i++)
- n_dimensions += (size_t) (dimensions_[i] > 1);
- return n_dimensions;
- }
-
-};
-
-} // namespace ISMRMRD
-#endif // ISMRMRD_NDARRAY_H
diff --git a/utilities/ismrmrd_phantom.cpp b/utilities/ismrmrd_phantom.cpp
index 2de7de8..dd9250c 100644
--- a/utilities/ismrmrd_phantom.cpp
+++ b/utilities/ismrmrd_phantom.cpp
@@ -12,11 +12,11 @@
namespace ISMRMRD {
-boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(std::vector<PhantomEllipse>& ellipses, unsigned int matrix_size)
+boost::shared_ptr<NDArray<complex_float_t> > phantom(std::vector<PhantomEllipse>& ellipses, unsigned int matrix_size)
{
- std::vector<unsigned int> dims(2,matrix_size);
- boost::shared_ptr<NDArrayContainer<std::complex<float> > > out(new NDArrayContainer< std::complex<float> >(dims));
- out->data_ = std::complex<float>(0.0,0.0);
+ std::vector<size_t> dims(2,matrix_size);
+ boost::shared_ptr<NDArray<complex_float_t> > out(new NDArray<complex_float_t>(dims));
+ memset(out->getData(), 0, out->getDataSize());
for (std::vector<PhantomEllipse>::iterator it = ellipses.begin(); it != ellipses.end(); it++) {
for (unsigned int y = 0; y < matrix_size; y++) {
float y_co = (1.0*y-(matrix_size>>1))/(matrix_size>>1);
@@ -24,7 +24,7 @@ boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(std::vector<
size_t index = y*matrix_size + x;
float x_co = (1.0*x-(matrix_size>>1))/(matrix_size>>1);
if (it->isInside(x_co,y_co)) {
- out->data_[index] += std::complex<float>(it->getAmplitude(),0.0);
+ out->getData()[index] += std::complex<float>(it->getAmplitude(),0.0);
}
}
}
@@ -33,7 +33,7 @@ boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(std::vector<
}
-boost::shared_ptr<NDArrayContainer< std::complex<float> > > shepp_logan_phantom(unsigned int matrix_size)
+boost::shared_ptr<NDArray<complex_float_t> > shepp_logan_phantom(unsigned int matrix_size)
{
boost::shared_ptr< std::vector<PhantomEllipse> > e = modified_shepp_logan_ellipses();
return phantom(*e, matrix_size);
@@ -73,15 +73,15 @@ boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses()
return out;
}
-boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius)
+boost::shared_ptr<NDArray<complex_float_t> > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius)
{
//This function is heavily inspired by the mri_birdcage.m Matlab script in Jeff Fessler's IRT packake
//http://web.eecs.umich.edu/~fessler/code/
- std::vector<unsigned int> dims(2,matrix_size);
+ std::vector<size_t> dims(2,matrix_size);
dims.push_back(ncoils);
- boost::shared_ptr<NDArrayContainer<std::complex<float> > > out(new NDArrayContainer< std::complex<float> >(dims));
- out->data_ = std::complex<float>(0.0,0.0);
+ boost::shared_ptr<NDArray<complex_float_t> > out(new NDArray<complex_float_t>(dims));
+ memset(out->getData(), 0, out->getDataSize());
for (int c = 0; c < ncoils; c++) {
float coilx = relative_radius*std::cos(c*(2*3.14159265359/ncoils));
@@ -94,7 +94,7 @@ boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_se
float x_co = (1.0*x-(matrix_size>>1))/(matrix_size>>1)-coilx;
float rr = std::sqrt(x_co*x_co+y_co*y_co);
float phi = atan2(x_co, -y_co) + coil_phase;
- out->data_[index] = std::polar(1 / rr, phi);
+ out->getData()[index] = std::polar(1 / rr, phi);
}
}
}
@@ -109,21 +109,21 @@ boost::mt19937& get_noise_seed()
return rng;
}
-int add_noise(NDArrayContainer< std::complex<float> >& a, float sd)
+int add_noise(NDArray<complex_float_t> & a, float sd)
{
boost::normal_distribution<float> nd(0.0, sd);
boost::variate_generator<boost::mt19937&,
boost::normal_distribution<float> > var_nor(get_noise_seed(), nd);
- for (size_t i = 0; i < a.data_.size(); i++) {
- a.data_[i] += std::complex<float>(var_nor(),var_nor());
+ for (size_t i = 0; i < a.getNumberOfElements(); i++) {
+ a.getData()[i] += std::complex<float>(var_nor(),var_nor());
}
return 0;
}
-int add_noise(ISMRMRD::Acquisition& a, float sd)
+int add_noise(Acquisition& a, float sd)
{
boost::normal_distribution<float> nd(0.0, sd);
@@ -137,7 +137,3 @@ int add_noise(ISMRMRD::Acquisition& a, float sd)
return 0;
}
};
-
-
-
-
diff --git a/utilities/ismrmrd_phantom.h b/utilities/ismrmrd_phantom.h
index da89b2e..fdae0ca 100644
--- a/utilities/ismrmrd_phantom.h
+++ b/utilities/ismrmrd_phantom.h
@@ -9,7 +9,6 @@
#include <vector>
#include <complex>
#include "ismrmrd/ismrmrd.h"
-#include "ismrmrd_ndarray.h"
#ifndef ISMRMRD_PHANTOM_H_
#define ISMRMRD_PHANTOM_H_
@@ -59,17 +58,12 @@ namespace ISMRMRD {
};
boost::shared_ptr< std::vector<PhantomEllipse> > shepp_logan_ellipses();
-
boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses();
-
- boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(NDArrayContainer<float>& coefficients, unsigned int matrix_size);
-
- boost::shared_ptr<NDArrayContainer< std::complex<float> > > shepp_logan_phantom(unsigned int matrix_size);
-
- boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius);
-
- int add_noise(NDArrayContainer< std::complex<float> >& a, float sd);
- int add_noise(ISMRMRD::Acquisition& a, float sd);
+ boost::shared_ptr<NDArray<complex_float_t> > phantom(std::vector<PhantomEllipse>& coefficients, unsigned int matrix_size);
+ boost::shared_ptr<NDArray<complex_float_t> > shepp_logan_phantom(unsigned int matrix_size);
+ boost::shared_ptr<NDArray<complex_float_t> > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius);
+ int add_noise(NDArray<complex_float_t> & a, float sd);
+ int add_noise(Acquisition & a, float sd);
};
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/ismrmrd.git
More information about the debian-science-commits
mailing list