[ismrmrd] 28/281: Complete first example of creating a dataset from scratch

Ghislain Vaillant ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:00:52 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 977a9365c9c003b905e6694742ed78dfc8cc37c5
Author: Michael Hansen <michael.hansen at nih.gov>
Date:   Wed Sep 5 21:37:26 2012 -0400

    Complete first example of creating a dataset from scratch
---
 CMakeLists.txt          |   3 +-
 test_create_dataset.cpp | 123 +++++++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 118 insertions(+), 8 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3eff982..ef124cb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -55,8 +55,9 @@ target_link_libraries(ismrmrd_test ${XERCESC_LIBRARIES} ismrmrd)
 find_package(FFTW3 COMPONENTS single)
 IF(FFTW3_FOUND)
 	MESSAGE("FFTW3 Found, building test applications")
+	INCLUDE_DIRECTORIES(${FFTW3_INDCLUDE_DIR})
 	add_executable(ismrmrd_create_dataset test_create_dataset.cpp)
-	target_link_libraries(ismrmrd_create_dataset ${XERCESC_LIBRARIES} ismrmrd)
+	target_link_libraries(ismrmrd_create_dataset ${XERCESC_LIBRARIES} ismrmrd ${FFTW3_LIBRARIES})
 ELSE(FFTW3_FOUND)
 	MESSAGE("FFTW3 NOT Found....cannot build test applications")
 ENDIF(FFTW3_FOUND)
diff --git a/test_create_dataset.cpp b/test_create_dataset.cpp
index aa051c0..3ceac5c 100644
--- a/test_create_dataset.cpp
+++ b/test_create_dataset.cpp
@@ -7,10 +7,11 @@
 
 #include <iostream>
 #include "ismrmrd_hdf5.h"
+#include "fftw3.h"
+#include "ismrmrd.hxx"
 
-using namespace ISMRMRD;
-
-template <typename T, int size_x, int size_y> int appendImageArray(IsmrmrdDataset& d, const char* varname)
+//Utility function for appending different sizes of arrays, used for testing here
+template <typename T, int size_x, int size_y> int appendImageArray(ISMRMRD::IsmrmrdDataset& d, const char* varname)
 {
 	T a[size_x*size_y];
 	std::vector<unsigned int> dims(2,0);
@@ -29,26 +30,42 @@ template <typename T, int size_x, int size_y> int appendImageArray(IsmrmrdDatase
 		}
 	}
 
-	NDArrayContainer<T> tmp(dims,a);
+	ISMRMRD::NDArrayContainer<T> tmp(dims,a);
 
 	return d.appendArray(tmp, varname);
 }
 
+//Helper function for the FFTW library
+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 main(int argc, char** argv)
 {
 	std::cout << "ISMRMRD Test Dataset Creation App" << std::endl;
 
-	IsmrmrdDataset d("testdata.h5","dataset");
+	const unsigned int readout = 256;
+	const unsigned int phase_encoding_lines = 128;
+
+	ISMRMRD::IsmrmrdDataset d("testdata.h5","dataset");
 
 	//Let's create the "original" image in the file for reference
-	if (appendImageArray< std::complex<float>, 256, 128 >(d, "the_square") < 0) {
+	if (appendImageArray< std::complex<float>, readout, phase_encoding_lines >(d, "the_square") < 0) {
 		std::cout << "Error adding image to dataset" << std::endl;
 		return -1;
 	}
 
 	//Read it back from the file
-	boost::shared_ptr< NDArrayContainer<std::complex<float> > > img_test =
+	boost::shared_ptr< ISMRMRD::NDArrayContainer<std::complex<float> > > img_test =
 			d.readArray< std::complex<float> >("the_square", 0);
 
 	if (img_test.get() == 0) {
@@ -62,6 +79,98 @@ int main(int argc, char** argv)
 	}
 	std::cout << std::endl;
 
+	//Let's FFT this image to k-space
+	fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*img_test->data_.size());
+
+	if (!tmp) {
+		std::cout << "Error allocating temporary storage for FFTW" << std::endl;
+		return -1;
+	}
+
+	fftshift(reinterpret_cast<std::complex<float>*>(tmp),&img_test->data_[0],img_test->dimensions_[0],img_test->dimensions_[1]);
+
+	//Create the FFTW plan
+	fftwf_plan p = fftwf_plan_dft_2d(img_test->dimensions_[0], img_test->dimensions_[1],
+	                                tmp,tmp, FFTW_FORWARD, FFTW_ESTIMATE);
+
+	fftwf_execute(p);
+
+	fftshift(&img_test->data_[0],reinterpret_cast<std::complex<float>*>(tmp),img_test->dimensions_[0],img_test->dimensions_[1]);
+
+	//Clean up.
+	fftwf_destroy_plan(p);
+	fftwf_free(tmp);
+
+	//Let keep the "original" k-space in the file for reference
+	if (d.appendArray(*img_test,"the_square_k") < 0) {
+		std::cout << "Error adding kspace to dataset" << std::endl;
+		return -1;
+	}
+
+	/* This is what one could to to get the magnitude of k-space
+	NDArrayContainer< float > abs_array;
+	abs_array.dimensions_ = img_test->dimensions_;
+	abs_array.data_.resize(img_test->data_.size(),0.0);
+
+	for (unsigned int i = 0; i < abs_array.data_.size();i++) {
+		abs_array.data_[i] = abs(img_test->data_[i]);
+	}
+
+	if (d.appendArray(abs_array,"the_square_k_abs") < 0) {
+		std::cout << "Error adding kspace to dataset" << std::endl;
+		return -1;
+	}
+	*/
+
+	//Let's append the data to the file
+	ISMRMRD::Acquisition acq;
+
+	acq.data_ = new float[readout*2];
+	if (!acq.data_) {
+		std::cout << "Error allocating memory for the acquisition" << std::endl;
+	}
+
+	for (unsigned int i = 0; i < phase_encoding_lines; i++) {
+		acq.head_.flags = 0;
+		//Set some flags
+		if (i == 0) {
+			acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_FIRST_IN_SLICE));
+		}
+		if (i == (phase_encoding_lines-1)) {
+			acq.setFlag(ISMRMRD::FlagBit(ISMRMRD::ACQ_LAST_IN_SLICE));
+		}
+		acq.head_.idx.kspace_encode_step_1 = i;
+		memcpy(acq.data_,&img_test->data_[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,phase_encoding_lines,1),ISMRMRD::fieldOfView_mm(600,300,6));
+	ISMRMRD::encodingSpaceType rs(ISMRMRD::matrixSize((readout>>1),phase_encoding_lines,1),ISMRMRD::fieldOfView_mm(300,300,6));
+	ISMRMRD::encodingLimitsType el;
+	el.kspace_encoding_step_1(ISMRMRD::limitType(0,phase_encoding_lines-1,(phase_encoding_lines>>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....
+
+	//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);
+
 	return 0;
 }
 

-- 
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