[ismrmrd] 121/281: Added first parts of utility functions for generating shepp-logan phantom datasets
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:01:05 UTC 2015
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to annotated tag ismrmrd0.5
in repository ismrmrd.
commit 1e49962f259e6eaf723f79cca79d0ddfcce8fa43
Author: Michael S. Hansen <michael.hansen at nih.gov>
Date: Mon Apr 1 22:54:35 2013 -0400
Added first parts of utility functions for generating shepp-logan phantom datasets
---
CMakeLists.txt | 2 +
utilities/CMakeLists.txt | 13 +++
utilities/generate_cartesian_shepp_logan.cpp | 126 +++++++++++++++++++++++++++
utilities/ismrmrd_fftw.h | 79 +++++++++++++++++
utilities/ismrmrd_phantom.cpp | 124 ++++++++++++++++++++++++++
utilities/ismrmrd_phantom.h | 73 ++++++++++++++++
utilities/ismrmrd_utilities_export.h | 21 +++++
7 files changed, 438 insertions(+)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 41bbfbc..9eeef74 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -69,3 +69,5 @@ INSTALL(FILES ${MatlabMFiles} DESTINATION matlab/+ismrmrd)
add_subdirectory(doc)
add_subdirectory(examples/c++)
add_subdirectory(tests)
+add_subdirectory(utilities)
+
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
new file mode 100644
index 0000000..74109d8
--- /dev/null
+++ b/utilities/CMakeLists.txt
@@ -0,0 +1,13 @@
+find_package(Boost REQUIRED)
+find_package(FFTW3 COMPONENTS single)
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/.. ${Boost_INCLUDE_DIR})
+
+ADD_LIBRARY(ismrmrd_utilities SHARED ismrmrd_phantom.cpp)
+INSTALL(TARGETS ismrmrd_utilities DESTINATION lib)
+
+IF(FFTW3_FOUND)
+ INCLUDE_DIRECTORIES(${FFTW3_INCLUDE_DIR})
+ ADD_EXECUTABLE(ismrmrd_generate_cartesian_shepp_logan generate_cartesian_shepp_logan.cpp)
+ TARGET_LINK_LIBRARIES(ismrmrd_generate_cartesian_shepp_logan ismrmrd_utilities ismrmrd ismrmrd_xsd ${FFTW3_LIBRARIES})
+ INSTALL(TARGETS ismrmrd_generate_cartesian_shepp_logan DESTINATION bin)
+ENDIF(FFTW3_FOUND)
\ No newline at end of file
diff --git a/utilities/generate_cartesian_shepp_logan.cpp b/utilities/generate_cartesian_shepp_logan.cpp
new file mode 100644
index 0000000..a094aa1
--- /dev/null
+++ b/utilities/generate_cartesian_shepp_logan.cpp
@@ -0,0 +1,126 @@
+/*
+ * generate_cartesian_shepp_logan.cpp
+ *
+ * Created on: Apr 1, 2013
+ * Author: Michael S. Hansen
+ *
+ */
+
+#include <iostream>
+#include "ismrmrd_phantom.h"
+#include "ismrmrd_hdf5.h"
+#include "ismrmrd_fftw.h"
+#include "ismrmrd.hxx"
+
+
+using namespace ISMRMRD;
+
+int main(int argc, char** argv)
+{
+ std::cout << "Generating Cartesian Shepp Logan Phantom" << std::endl;
+
+ unsigned int matrix_size = 256; //Matrix size
+ unsigned int ncoils = 8; //Number of coils
+ unsigned int ros = 2; //Readout ovesampling
+ unsigned int repetitions = 10;
+ float noise_level = 0.05;
+
+ 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);
+
+ std::vector<unsigned int> 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);
+
+ for (unsigned int c = 0; c < ncoils; c++) {
+ for (unsigned int y = 0; y < matrix_size; y++) {
+ for (unsigned int x = 0; x < matrix_size; x++) {
+ 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];
+ }
+ }
+ }
+
+
+ //Let's append the data to the file
+ IsmrmrdDataset d("testdata.h5","dataset");
+ Acquisition acq;
+ acq.setActiveChannels(ncoils);
+ acq.setAvailableChannels(ncoils);
+ size_t readout = matrix_size*ros;
+ acq.setNumberOfSamples(readout);
+ acq.setCenterSample(readout>>1);
+
+ for (unsigned int r = 0; r < repetitions; r++) {
+ NDArrayContainer<std::complex<float> > cm = coil_images;
+ fft2c(cm);
+ add_noise(cm,noise_level);
+ for (unsigned int i = 0; i < matrix_size; i++) {
+ acq.setFlags(0);
+
+ //Set some flags
+ if (i == 0) {
+ acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_FIRST_IN_SLICE));
+ }
+ if (i == (matrix_size-1)) {
+ acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_LAST_IN_SLICE));
+ }
+ acq.getIdx().kspace_encode_step_1 = i;
+ acq.getIdx().repetition = r;
+ acq.setSampleTimeUs(5.0);
+ for (unsigned int c = 0; c < ncoils; c++) {
+ memcpy(&acq[c*readout*2],&(cm.data_[c*matrix_size*readout + i*readout]),sizeof(float)*readout*2);
+ }
+ d.appendAcquisition(&acq);
+ }
+ }
+
+ //Let's create a header, we will use the C++ class generated by XSD
+ ISMRMRD::experimentalConditionsType exp(63500000); //~1.5T
+ ISMRMRD::ismrmrdHeader h(exp);
+
+ //Create an encoding section
+ ISMRMRD::encodingSpaceType es(ISMRMRD::matrixSize(readout,matrix_size,1),ISMRMRD::fieldOfView_mm(600,300,6));
+ ISMRMRD::encodingSpaceType rs(ISMRMRD::matrixSize((readout>>1),matrix_size,1),ISMRMRD::fieldOfView_mm(300,300,6));
+ ISMRMRD::encodingLimitsType el;
+ el.kspace_encoding_step_1(ISMRMRD::limitType(0,matrix_size-1,(matrix_size>>1)));
+ ISMRMRD::encoding e(es,rs,el,ISMRMRD::trajectoryType::cartesian);
+
+ //Add the encoding section to the header
+ h.encoding().push_back(e);
+
+ //Add any additional fields that you may want would go here....
+
+ //e.g. parallel imaging
+ //ISMRMRD::parallelImagingType parallel(ISMRMRD::accelerationFactorType(2,1));
+ //parallel.calibrationMode(ISMRMRD::calibrationModeType::embedded);
+ //h.parallelImaging(parallel);
+
+ //Serialize the header
+ xml_schema::namespace_infomap map;
+ map[""].name = "http://www.ismrm.org/ISMRMRD";
+ map[""].schema = "ismrmrd.xsd";
+ std::stringstream str;
+ ISMRMRD::ismrmrdHeader_(str, h, map);
+ std::string xml_header = str.str();
+
+ //Write the header to the data file.
+ d.writeHeader(xml_header);
+
+
+ d.appendArray(*phantom,"phantom");
+ d.appendArray(*coils,"csm");
+ d.appendArray(coil_images,"coil_images");
+
+ return 0;
+}
+
+
+
+
diff --git a/utilities/ismrmrd_fftw.h b/utilities/ismrmrd_fftw.h
new file mode 100644
index 0000000..c67effc
--- /dev/null
+++ b/utilities/ismrmrd_fftw.h
@@ -0,0 +1,79 @@
+/*
+ * ismrmrd_fftw.h
+ *
+ * Created on: Apr 1, 2013
+ * Author: Michael S. Hansen
+ */
+
+#include "fftw3.h"
+
+namespace ISMRMRD {
+
+template<typename TI, typename TO> void circshift(TO *out, const TI *in, int xdim, int ydim, int xshift, int yshift)
+{
+ for (int i =0; i < ydim; i++) {
+ int ii = (i + yshift) % ydim;
+ for (int j = 0; j < xdim; j++) {
+ int jj = (j + xshift) % xdim;
+ out[ii * xdim + jj] = in[i * xdim + j];
+ }
+ }
+}
+
+#define fftshift(out, in, x, y) circshift(out, in, x, y, (x/2), (y/2))
+
+
+int fft2c(NDArrayContainer<std::complex<float> >& a, bool forward)
+{
+ if (a.dimensions_.size() < 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;
+
+ //Array for transformation
+ fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*a.data_.size());
+
+ 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]);
+
+ //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);
+ } else {
+ p = fftwf_plan_dft_2d(a.dimensions_[1], a.dimensions_[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]);
+
+ //Clean up.
+ fftwf_destroy_plan(p);
+ }
+
+ std::complex<float> scale(std::sqrt(1.0f*elements),0.0);
+ a.data_ /= scale;
+
+ fftwf_free(tmp);
+ return 0;
+}
+
+int fft2c(NDArrayContainer<std::complex<float> >& a) {
+ return fft2c(a,true);
+}
+
+int ifft2c(NDArrayContainer<std::complex<float> >& a) {
+ return fft2c(a,false);
+}
+
+};
diff --git a/utilities/ismrmrd_phantom.cpp b/utilities/ismrmrd_phantom.cpp
new file mode 100644
index 0000000..5472f83
--- /dev/null
+++ b/utilities/ismrmrd_phantom.cpp
@@ -0,0 +1,124 @@
+/*
+ * ismrmrd_phanthom.cpp
+ *
+ * Created on: Apr 1, 2013
+ * Author: Michael S. Hansen
+ */
+
+#include "ismrmrd_phantom.h"
+#include <boost/random.hpp>
+#include <boost/random/normal_distribution.hpp>
+
+namespace ISMRMRD {
+
+
+boost::shared_ptr<NDArrayContainer< std::complex<float> > > 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);
+ 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);
+ for (unsigned int x = 0; x < matrix_size; x++) {
+ 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);
+ }
+ }
+ }
+ }
+ return out;
+}
+
+
+boost::shared_ptr<NDArrayContainer< std::complex<float> > > shepp_logan_phantom(unsigned int matrix_size)
+{
+ boost::shared_ptr< std::vector<PhantomEllipse> > e = modified_shepp_logan_ellipses();
+ return phantom(*e, matrix_size);
+}
+
+
+boost::shared_ptr< std::vector<PhantomEllipse> > shepp_logan_ellipses()
+{
+ boost::shared_ptr< std::vector<PhantomEllipse> > out(new std::vector<PhantomEllipse>);
+ out->push_back(PhantomEllipse(1, .69, .92, 0, 0, 0));
+ out->push_back(PhantomEllipse(-.98, .6624, .8740, 0, -.0184, 0));
+ out->push_back(PhantomEllipse(-.02, .1100, .3100, .22, 0, -18));
+ out->push_back(PhantomEllipse(-.02, .1600, .4100, -.22, 0, 18));
+ out->push_back(PhantomEllipse(.01, .2100, .2500, 0, .35, 0));
+ out->push_back(PhantomEllipse(.01, .0460, .0460, 0, .1, 0));
+ out->push_back(PhantomEllipse(.01, .0460, .0460, 0, -.1, 0));
+ out->push_back(PhantomEllipse(.01, .0460, .0230, -.08, -.605, 0));
+ out->push_back(PhantomEllipse(.01, .0230, .0230, 0, -.606, 0));
+ out->push_back(PhantomEllipse( .01, .0230, .0460, .06, -.605, 0));
+
+ return out;
+}
+
+boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses()
+{
+ boost::shared_ptr< std::vector<PhantomEllipse> > out(new std::vector<PhantomEllipse>);
+ out->push_back(PhantomEllipse( 1, .69, .92, 0, 0, 0));
+ out->push_back(PhantomEllipse(-.8, .6624, .8740, 0, -.0184, 0));
+ out->push_back(PhantomEllipse(-.2, .1100, .3100, .22, 0, -18));
+ out->push_back(PhantomEllipse(-.2, .1600, .4100, -.22, 0, 18));
+ out->push_back(PhantomEllipse(.1, .2100, .2500, 0, .35, 0));
+ out->push_back(PhantomEllipse(.1, .0460, .0460, 0, .1, 0));
+ out->push_back(PhantomEllipse(.1, .0460, .0460, 0, -.1, 0));
+ out->push_back(PhantomEllipse(.1, .0460, .0230, -.08, -.605, 0));
+ out->push_back(PhantomEllipse(.1, .0230, .0230, 0, -.606, 0));
+ out->push_back(PhantomEllipse( .1, .0230, .0460, .06, -.605, 0 ));
+ return out;
+}
+
+boost::shared_ptr<NDArrayContainer< std::complex<float> > > 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);
+ 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);
+
+ for (int c = 0; c < ncoils; c++) {
+ float coilx = relative_radius*std::cos(c*(2*3.14159265359/ncoils));
+ float coily = relative_radius*std::sin(c*(2*3.14159265359/ncoils));
+ float coil_phase = -c*(2*3.14159265359/ncoils);
+ for (unsigned int y = 0; y < matrix_size; y++) {
+ float y_co = (1.0*y-(matrix_size>>1))/(matrix_size>>1)-coily;
+ for (unsigned int x = 0; x < matrix_size; x++) {
+ size_t index = c*matrix_size*matrix_size+ y*matrix_size + x;
+ 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);
+ }
+ }
+ }
+
+ return out;
+}
+
+int add_noise(NDArrayContainer< std::complex<float> >& a, float sd)
+{
+
+ static boost::mt19937 rng;
+ boost::normal_distribution<float> nd(0.0, sd);
+ boost::variate_generator<boost::mt19937&,
+ boost::normal_distribution<float> > var_nor(rng, nd);
+
+ for (size_t i = 0; i < a.data_.size(); i++) {
+ a.data_[i] += std::complex<float>(var_nor(),var_nor());
+ }
+
+ return 0;
+}
+
+};
+
+
+
+
diff --git a/utilities/ismrmrd_phantom.h b/utilities/ismrmrd_phantom.h
new file mode 100644
index 0000000..3a4a827
--- /dev/null
+++ b/utilities/ismrmrd_phantom.h
@@ -0,0 +1,73 @@
+/*
+ * ismrmrd_phantom.h
+ *
+ * Created on: Apr 1, 2013
+ * Author: Michael S. Hansen
+ */
+
+#include <boost/shared_ptr.hpp>
+#include <complex>
+#include "ismrmrd.h"
+#include "ismrmrd_utilities_export.h"
+
+#ifndef ISMRMRD_PHANTOM_H_
+#define ISMRMRD_PHANTOM_H_
+
+namespace ISMRMRD {
+
+ class PhantomEllipse
+ {
+ public:
+ PhantomEllipse(float A, float a, float b, float x0, float y0, float phi)
+ : A_(A)
+ , a_(a)
+ , b_(b)
+ , x0_(x0)
+ , y0_(y0)
+ , phi_(phi)
+ {
+
+ }
+
+ bool isInside(float x, float y)
+ {
+ float asq = a_*a_; // a^2
+ float bsq = b_*b_; // b^2
+ float phi = phi_*3.14159265359/180; // rotation angle in radians
+ float x0 = x-x0_; // x offset
+ float y0 = y-y0_; // y offset
+ float cosp = cos(phi);
+ float sinp = sin(phi);
+ return (((x0*cosp + y0*sinp)*(x0*cosp + y0*sinp))/asq + ((y0*cosp - x0*sinp)*(y0*cosp - x0*sinp))/bsq <= 1);
+ }
+
+ float getAmplitude()
+ {
+ return A_;
+ }
+
+
+ protected:
+
+ float A_;
+ float a_;
+ float b_;
+ float x0_;
+ float y0_;
+ float phi_;
+ };
+
+ EXPORTISMRMRDUTILITIES boost::shared_ptr< std::vector<PhantomEllipse> > shepp_logan_ellipses();
+ EXPORTISMRMRDUTILITIES boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses();
+
+ EXPORTISMRMRDUTILITIES boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(NDArrayContainer<float>& coefficients, unsigned int matrix_size);
+ EXPORTISMRMRDUTILITIES boost::shared_ptr<NDArrayContainer< std::complex<float> > > shepp_logan_phantom(unsigned int matrix_size);
+
+ EXPORTISMRMRDUTILITIES boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius);
+
+ EXPORTISMRMRDUTILITIES int add_noise(NDArrayContainer< std::complex<float> >& a, float sd);
+
+
+};
+
+#endif /* ISMRMRD_PHANTOM_H_ */
diff --git a/utilities/ismrmrd_utilities_export.h b/utilities/ismrmrd_utilities_export.h
new file mode 100644
index 0000000..e9101a0
--- /dev/null
+++ b/utilities/ismrmrd_utilities_export.h
@@ -0,0 +1,21 @@
+/*
+ * ismrmrd_utilities_export.h
+ *
+ * Created on: Apr 1, 2013
+ * Author: Michael S. Hansen
+ */
+
+#ifndef ISMRMRD_UTILITIES_EXPORT_H_
+#define ISMRMRD_UTILITIES_EXPORT_H_
+
+#if defined (WIN32)
+#if defined (ismrmrd_utilities_EXPORTS)
+#define EXPORTISMRMRDUTILITIES __declspec(dllexport)
+#else
+#define EXPORTISMRMRDUTILITIES __declspec(dllimport)
+#endif
+#else
+#define EXPORTISMRMRDUTILITIES
+#endif
+
+#endif /* ISMRMRD_UTILITIES_EXPORT_H_ */
--
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