[ismrmrd] 122/281: Further improvements to shepp-logan phantom generator
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 643bf27b6ec1eda508ff64adb576a602ee072928
Author: Michael S. Hansen <michael.hansen at nih.gov>
Date: Tue Apr 2 09:46:30 2013 -0400
Further improvements to shepp-logan phantom generator
---
utilities/CMakeLists.txt | 4 +-
utilities/generate_cartesian_shepp_logan.cpp | 139 +++++++++++++++++++++------
utilities/ismrmrd_phantom.cpp | 23 ++++-
utilities/ismrmrd_phantom.h | 2 +-
4 files changed, 134 insertions(+), 34 deletions(-)
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index 74109d8..d29b62c 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -1,4 +1,4 @@
-find_package(Boost REQUIRED)
+find_package(Boost COMPONENTS program_options REQUIRED REQUIRED)
find_package(FFTW3 COMPONENTS single)
INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/.. ${Boost_INCLUDE_DIR})
@@ -8,6 +8,6 @@ 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})
+ TARGET_LINK_LIBRARIES(ismrmrd_generate_cartesian_shepp_logan ismrmrd_utilities ismrmrd ismrmrd_xsd ${FFTW3_LIBRARIES} ${Boost_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
index a094aa1..f875330 100644
--- a/utilities/generate_cartesian_shepp_logan.cpp
+++ b/utilities/generate_cartesian_shepp_logan.cpp
@@ -12,18 +12,58 @@
#include "ismrmrd_fftw.h"
#include "ismrmrd.hxx"
+#include <boost/program_options.hpp>
using namespace ISMRMRD;
+namespace po = boost::program_options;
int main(int argc, char** argv)
{
- std::cout << "Generating Cartesian Shepp Logan Phantom" << std::endl;
+ /** TODO
+ *
+ * Noise samples
+ * Acceleration
+ * k-space coordinates
+ *
+ */
+
+ unsigned int matrix_size; //Matrix size
+ unsigned int ncoils; //Number of coils
+ unsigned int ros; //Readout ovesampling
+ unsigned int repetitions;
+ unsigned int acc_factor;
+ float noise_level;
+ std::string outfile;
+ std::string dataset;
+ bool store_coordinates;
+ bool noise_calibration;
+
+ po::options_description desc("Allowed options");
+ desc.add_options()
+ ("help,h", "produce help message")
+ ("matrix,m", po::value<unsigned int>(&matrix_size)->default_value(256), "Matrix Size")
+ ("coils,c", po::value<unsigned int>(&ncoils)->default_value(8), "Number of Coils")
+ ("oversampling,O", po::value<unsigned int>(&ros)->default_value(2), "Readout oversampling")
+ ("repetitions,r", po::value<unsigned int>(&repetitions)->default_value(1), "Repetitions")
+ ("acceleration,a", po::value<unsigned int>(&acc_factor)->default_value(1), "Acceleration factor")
+ ("noise-level,n", po::value<float>(&noise_level)->default_value(0.05f,"0.05"), "Noise Level")
+ ("output,o", po::value<std::string>(&outfile)->default_value("testdata.h5"), "Output File Name")
+ ("dataset,d", po::value<std::string>(&dataset)->default_value("dataset"), "Output Dataset Name")
+ ("noise-calibration,C", po::value<bool>(&noise_calibration)->zero_tokens(), "Add noise calibration")
+ ("k-coordinates,k", po::value<bool>(&store_coordinates)->zero_tokens(), "Store k-space coordinates")
+ ;
+
+ po::variables_map vm;
+ po::store(po::parse_command_line(argc, argv, desc), vm);
+ po::notify(vm);
+
+ if (vm.count("help")) {
+ std::cout << desc << "\n";
+ return 1;
+ }
- 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;
+ std::cout << "Generating Cartesian Shepp Logan Phantom" << std::endl;
+ std::cout << "Accelleration: " << 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);
@@ -49,7 +89,7 @@ int main(int argc, char** argv)
//Let's append the data to the file
- IsmrmrdDataset d("testdata.h5","dataset");
+ IsmrmrdDataset d(outfile.c_str(),dataset.c_str());
Acquisition acq;
acq.setActiveChannels(ncoils);
acq.setAvailableChannels(ncoils);
@@ -57,39 +97,78 @@ int main(int argc, char** argv)
acq.setNumberOfSamples(readout);
acq.setCenterSample(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);
+ }
+
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));
+ 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);
}
- 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::acquisitionSystemInformationType sys;
+
+ sys.institutionName("ISMRM Synthetic Imaging Lab");
+ sys.receiverChannels(ncoils);
+
+
ISMRMRD::ismrmrdHeader h(exp);
+ 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);
//Add the encoding section to the header
@@ -98,9 +177,11 @@ int main(int argc, char** argv)
//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);
+ if (acc_factor > 1) {
+ ISMRMRD::parallelImagingType parallel(ISMRMRD::accelerationFactorType(acc_factor,1));
+ parallel.calibrationMode(ISMRMRD::calibrationModeType::interleaved);
+ h.parallelImaging(parallel);
+ }
//Serialize the header
xml_schema::namespace_infomap map;
diff --git a/utilities/ismrmrd_phantom.cpp b/utilities/ismrmrd_phantom.cpp
index 5472f83..9959254 100644
--- a/utilities/ismrmrd_phantom.cpp
+++ b/utilities/ismrmrd_phantom.cpp
@@ -102,13 +102,19 @@ boost::shared_ptr<NDArrayContainer< std::complex<float> > > generate_birdcage_se
return out;
}
+
+boost::mt19937& get_noise_seed()
+{
+ static boost::mt19937 rng;
+ return rng;
+}
+
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);
+ 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());
@@ -117,6 +123,19 @@ int add_noise(NDArrayContainer< std::complex<float> >& a, float sd)
return 0;
}
+int add_noise(ISMRMRD::Acquisition& 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.getData().size(); i++) {
+ a[i] += var_nor();
+ }
+
+ return 0;
+}
};
diff --git a/utilities/ismrmrd_phantom.h b/utilities/ismrmrd_phantom.h
index 3a4a827..fc0b65b 100644
--- a/utilities/ismrmrd_phantom.h
+++ b/utilities/ismrmrd_phantom.h
@@ -66,7 +66,7 @@ namespace ISMRMRD {
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);
-
+ EXPORTISMRMRDUTILITIES int add_noise(ISMRMRD::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