[ismrmrd] 85/177: Working on shepp-logan phantom utility.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:02:05 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 c0c842a9e8315166c719ab53bdcab497e808584c
Author: Souheil Inati <souheil.inati at nih.gov>
Date: Tue Sep 23 20:55:54 2014 -0400
Working on shepp-logan phantom utility.
---
CMakeLists.txt | 6 +-
utilities/CMakeLists.txt | 39 ++---
utilities/generate_cartesian_shepp_logan.cpp | 179 +++++++++++----------
utilities/ismrmrd_fftw.h | 1 +
utilities/ismrmrd_ndarray.h | 91 +++++++++++
utilities/ismrmrd_phantom.cpp | 4 +-
utilities/ismrmrd_phantom.h | 20 +--
utilities/ismrmrd_utilities_export.h | 21 ---
{examples/c++ => utilities}/test_recon_dataset.cpp | 0
9 files changed, 220 insertions(+), 141 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b02ad41..7b6a6c5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
project(ISMRMRD)
# add project specific cmake find modules
-list(APPEND CMAKE_MODULE_PATH cmake)
+list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake)
# set the build type to Release if not specified
IF(NOT CMAKE_BUILD_TYPE)
@@ -133,9 +133,9 @@ enable_testing()
# process subdirectories
#add_subdirectory(examples/c++)
#add_subdirectory(examples/c)
-#add_subdirectory(utilities)
+add_subdirectory(utilities)
#add_subdirectory(tests)
#add_subdirectory(doc)
#add_subdirectory(matlab)
#add_subdirectory(bindings)
-add_subdirectory(info)
\ No newline at end of file
+add_subdirectory(info)
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index 003adbd..d7f0560 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -2,28 +2,23 @@ find_package(Boost COMPONENTS program_options)
find_package(FFTW3 COMPONENTS single)
if(FFTW3_FOUND)
- message("FFTW3 Found, building utilities")
- include_directories(${ISMRMRD_SOURCE_DIR} ${ISMRMRD_SCHEMA_SOURCE_DIR}
- ${Boost_INCLUDE_DIR} ${FFTW3_INCLUDE_DIR})
- add_library(ismrmrd_utilities SHARED ismrmrd_phantom.cpp)
- install(TARGETS ismrmrd_utilities DESTINATION ${ISMRMRD_INSTALL_LIB_DIR})
-
- if(WIN32)
- set_source_files_properties(${XSDS_SOURCES} PROPERTIES GENERATED TRUE)
- add_executable(ismrmrd_generate_cartesian_shepp_logan
- generate_cartesian_shepp_logan.cpp ${XSDS_SOURCES})
- target_link_libraries(ismrmrd_generate_cartesian_shepp_logan
- ismrmrd ismrmrd_utilities ${Boost_LIBRARIES} ${FFTW3_LIBRARIES}
- ${XERCESC_LIBRARIES})
- else(WIN32)
- add_executable(ismrmrd_generate_cartesian_shepp_logan
- generate_cartesian_shepp_logan.cpp)
- target_link_libraries(ismrmrd_generate_cartesian_shepp_logan
- ismrmrd ismrmrd_utilities ismrmrd_xsd
- ${Boost_PROGRAM_OPTIONS_LIBRARY} ${FFTW3_LIBRARIES})
- endif(WIN32)
- install(TARGETS ismrmrd_generate_cartesian_shepp_logan
- DESTINATION ${ISMRMRD_INSTALL_BIN_DIR})
+ message("FFTW3 Found, building utilities")
+ include_directories(
+ ${CMAKE_SOURCE_DIR/include}
+ ${Boost_INCLUDE_DIR}
+ ${FFTW3_INCLUDE_DIR})
+
+ # Shepp-Logan phantom
+ add_executable(ismrmrd_generate_cartesian_shepp_logan
+ generate_cartesian_shepp_logan.cpp
+ ismrmrd_phantom.cpp)
+ target_link_libraries( ismrmrd_generate_cartesian_shepp_logan
+ ismrmrd
+ ${Boost_PROGRAM_OPTIONS_LIBRARY}
+ ${FFTW3_LIBRARIES})
+ install(TARGETS ismrmrd_generate_cartesian_shepp_logan
+ DESTINATION bin)
+
else(FFTW3_FOUND)
message("FFTW3 NOT Found, cannot build utilities")
endif(FFTW3_FOUND)
diff --git a/utilities/generate_cartesian_shepp_logan.cpp b/utilities/generate_cartesian_shepp_logan.cpp
index d994441..98e2fec 100644
--- a/utilities/generate_cartesian_shepp_logan.cpp
+++ b/utilities/generate_cartesian_shepp_logan.cpp
@@ -7,10 +7,11 @@
*/
#include <iostream>
+#include "ismrmrd/ismrmrd.h"
+#include "ismrmrd/xml.h"
+#include "ismrmrd/dataset.h"
#include "ismrmrd_phantom.h"
-#include "ismrmrd_hdf5.h"
#include "ismrmrd_fftw.h"
-#include "ismrmrd.hxx"
#include <boost/program_options.hpp>
@@ -89,88 +90,99 @@ int main(int argc, char** argv)
//Let's append the data to the file
- IsmrmrdDataset d(outfile.c_str(),dataset.c_str());
+ //Create if needed
+ Dataset d(outfile.c_str(),dataset.c_str(), true);
Acquisition acq;
- acq.setActiveChannels(ncoils);
- acq.setAvailableChannels(ncoils);
+ acq.available_channels() = ncoils;
+ acq.active_channels(ncoils);
size_t readout = matrix_size*ros;
- acq.setNumberOfSamples(readout);
- acq.setCenterSample(readout>>1);
+ acq.number_of_samples(readout);
+ acq.center_sample() = (readout>>1);
if (noise_calibration) {
- acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_IS_NOISE_MEASUREMENT));
- acq.setNumberOfSamples(readout);
- for (size_t s = 0; s < acq.getData().size();s++) {
- acq[s] = 0.0;
- }
- add_noise(acq,noise_level);
- acq.setSampleTimeUs(5.0);
- d.appendAcquisition(&acq);
- }
-
- std::valarray<float> k;
- if (store_coordinates) {
- k.resize(2*readout);
+ acq.setFlag(ISMRMRD_ACQ_IS_NOISE_MEASUREMENT);
+ acq.number_of_samples(readout);
+ for (uint64_t s = 0; s < acq.number_of_samples()*acq.active_channels(); s++) {
+ acq.getData()[s] = 0.0;
+ }
+ add_noise(acq,noise_level);
+ acq.sample_time_us() = 5.0;
+ d.appendAcquisition(acq);
}
for (unsigned int r = 0; r < repetitions; r++) {
- for (unsigned int a = 0; a < acc_factor; a++) {
- NDArrayContainer<std::complex<float> > cm = coil_images;
- fft2c(cm);
- add_noise(cm,noise_level);
- for (int i = a; i < matrix_size; i+=acc_factor) {
- acq.setFlags(0);
- acq.setNumberOfSamples(readout);
-
- //Set some flags
- if (i == a) {
- acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_FIRST_IN_SLICE));
- }
- if (i >= (matrix_size-acc_factor)) {
- acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_LAST_IN_SLICE));
- }
- acq.getIdx().kspace_encode_step_1 = i;
- acq.getIdx().repetition = r*acc_factor + a;
- 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);
- }
-
- if (store_coordinates) {
- acq.setTrajectoryDimensions(2);
- float ky = (1.0*i-(matrix_size>>1))/(1.0*matrix_size);
- for (int x = 0; x < readout; x++) {
- float kx = (1.0*x-(readout>>1))/(1.0*readout);
- k[x*2 ] = kx;
- k[x*2+1] = ky;
- }
- acq.setTraj(k);
- }
- d.appendAcquisition(&acq);
- }
- }
+ for (unsigned int a = 0; a < acc_factor; a++) {
+ NDArrayContainer<std::complex<float> > cm = coil_images;
+ fft2c(cm);
+ add_noise(cm,noise_level);
+ for (int i = a; i < matrix_size; i+=acc_factor) {
+ acq.clearAllFlags();
+ acq.number_of_samples(readout);
+
+ //Set some flags
+ if (i == a) {
+ acq.setFlag(ISMRMRD_ACQ_FIRST_IN_SLICE);
+ }
+ if (i >= (matrix_size-acc_factor)) {
+ acq.setFlag(ISMRMRD_ACQ_LAST_IN_SLICE);
+ }
+ acq.idx().kspace_encode_step_1 = i;
+ acq.idx().repetition = r*acc_factor + a;
+ 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]),
+ sizeof(complex_float_t)*readout);
+ }
+
+ if (store_coordinates) {
+ acq.trajectory_dimensions(2);
+ float ky = (1.0*i-(matrix_size>>1))/(1.0*matrix_size);
+ for (int x = 0; x < readout; x++) {
+ float kx = (1.0*x-(readout>>1))/(1.0*readout);
+ acq.getTraj()[x*2 ] = kx;
+ acq.getTraj()[x*2+1] = ky;
+ }
+ }
+ d.appendAcquisition(acq);
+ }
+ }
}
- //Let's create a header, we will use the C++ class generated by XSD
- ISMRMRD::experimentalConditionsType exp(63500000); //~1.5T
- ISMRMRD::acquisitionSystemInformationType sys;
-
- sys.institutionName("ISMRM Synthetic Imaging Lab");
- sys.receiverChannels(ncoils);
+ //Let's create a header, we will use the C++ classes in ismrmrd/xml.h
+ IsmrmrdHeader h;
+ h.experimentalConditions.H1resonanceFrequency_Hz = 63500000; //~1.5T
-
- ISMRMRD::ismrmrdHeader h(exp);
-
- h.acquisitionSystemInformation(sys);
+ AcquisitionSystemInformation sys;
+ sys.institutionName = "ISMRM Synthetic Imaging Lab";
+ sys.receiverChannels = ncoils;
+ h.acquisitionSystemInformation = sys;
//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)));
- el.repetition(ISMRMRD::limitType(0,repetitions*acc_factor,0));
- ISMRMRD::encoding e(es,rs,el,ISMRMRD::trajectoryType::cartesian);
-
+ Encoding e;
+ e.encodedSpace.matrixSize.x = readout;
+ e.encodedSpace.matrixSize.y = matrix_size;
+ e.encodedSpace.matrixSize.z = 1;
+ e.encodedSpace.fieldOfView_mm.x = 600;
+ e.encodedSpace.fieldOfView_mm.y = 300;
+ e.encodedSpace.fieldOfView_mm.z = 6;
+ e.reconSpace.matrixSize.x = readout/2;
+ e.reconSpace.matrixSize.y = matrix_size;
+ e.reconSpace.matrixSize.z = 1;
+ e.reconSpace.fieldOfView_mm.x = 300;
+ e.reconSpace.fieldOfView_mm.y = 300;
+ e.reconSpace.fieldOfView_mm.z = 6;
+ TrajectoryDescription td;
+ td.identifier = "cartesian";
+ e.trajectoryDescription = td;
+
+ //encodingLimitsType el;
+ //el.kspace_encoding_step_1(limitType(0,matrix_size-1,(matrix_size>>1)));
+ //el.repetition(limitType(0,repetitions*acc_factor,0));
+ //encoding e(es,rs,el,trajectoryType::cartesian);
+
+#ifdef YOMAMA
+
//e.g. parallel imaging
if (acc_factor > 1) {
ISMRMRD::parallelImagingType parallel(ISMRMRD::accelerationFactorType(acc_factor,1));
@@ -178,27 +190,26 @@ int main(int argc, char** argv)
e.parallelImaging(parallel);
}
+#endif // YOMAMA
+
//Add the encoding section to the header
- h.encoding().push_back(e);
+ h.encoding.push_back(e);
//Add any additional fields that you may want would go here....
-
//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();
-
+ std::stringstream str;
+ ISMRMRD::serialize( h, str);
+ std::string xml_header = str.str();
+ std::cout << xml_header << std::endl;
+
//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");
+ //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
index c67effc..94e68b1 100644
--- a/utilities/ismrmrd_fftw.h
+++ b/utilities/ismrmrd_fftw.h
@@ -6,6 +6,7 @@
*/
#include "fftw3.h"
+#include "ismrmrd_ndarray.h"
namespace ISMRMRD {
diff --git a/utilities/ismrmrd_ndarray.h b/utilities/ismrmrd_ndarray.h
new file mode 100644
index 0000000..f37c57f
--- /dev/null
+++ b/utilities/ismrmrd_ndarray.h
@@ -0,0 +1,91 @@
+/**
+ * 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 9959254..2de7de8 100644
--- a/utilities/ismrmrd_phantom.cpp
+++ b/utilities/ismrmrd_phantom.cpp
@@ -130,8 +130,8 @@ int add_noise(ISMRMRD::Acquisition& a, float sd)
boost::variate_generator<boost::mt19937&,
boost::normal_distribution<float> > var_nor(get_noise_seed(), nd);
- for (size_t i = 0; i < a.getData().size(); i++) {
- a[i] += var_nor();
+ for (size_t i = 0; i < a.number_of_samples()*a.active_channels(); i++) {
+ a.getData()[i] += std::complex<float>(var_nor(), var_nor());
}
return 0;
diff --git a/utilities/ismrmrd_phantom.h b/utilities/ismrmrd_phantom.h
index 72362aa..da89b2e 100644
--- a/utilities/ismrmrd_phantom.h
+++ b/utilities/ismrmrd_phantom.h
@@ -8,8 +8,8 @@
#include <boost/shared_ptr.hpp>
#include <vector>
#include <complex>
-#include "ismrmrd.h"
-#include "ismrmrd_utilities_export.h"
+#include "ismrmrd/ismrmrd.h"
+#include "ismrmrd_ndarray.h"
#ifndef ISMRMRD_PHANTOM_H_
#define ISMRMRD_PHANTOM_H_
@@ -58,16 +58,18 @@ namespace ISMRMRD {
float phi_;
};
- EXPORTISMRMRDUTILITIES boost::shared_ptr< std::vector<PhantomEllipse> > shepp_logan_ellipses();
- EXPORTISMRMRDUTILITIES boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses();
+ boost::shared_ptr< std::vector<PhantomEllipse> > 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);
+ boost::shared_ptr< std::vector<PhantomEllipse> > modified_shepp_logan_ellipses();
- EXPORTISMRMRDUTILITIES boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_sensititivies(unsigned int matrix_size, unsigned int ncoils, float relative_radius);
+ boost::shared_ptr<NDArrayContainer< std::complex<float> > > phantom(NDArrayContainer<float>& coefficients, unsigned int matrix_size);
- EXPORTISMRMRDUTILITIES int add_noise(NDArrayContainer< std::complex<float> >& a, float sd);
- EXPORTISMRMRDUTILITIES int add_noise(ISMRMRD::Acquisition& a, float sd);
+ 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);
};
diff --git a/utilities/ismrmrd_utilities_export.h b/utilities/ismrmrd_utilities_export.h
deleted file mode 100644
index e9101a0..0000000
--- a/utilities/ismrmrd_utilities_export.h
+++ /dev/null
@@ -1,21 +0,0 @@
-/*
- * 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_ */
diff --git a/examples/c++/test_recon_dataset.cpp b/utilities/test_recon_dataset.cpp
similarity index 100%
rename from examples/c++/test_recon_dataset.cpp
rename to utilities/test_recon_dataset.cpp
--
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