[ismrmrd] 110/177: Fixed recon2d example utility.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:02:08 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 c4c5ac9170cba11fe98c095bc4b41545b133e66a
Author: Souheil Inati <souheil.inati at nih.gov>
Date: Sat Sep 27 00:02:38 2014 -0400
Fixed recon2d example utility.
---
utilities/CMakeLists.txt | 16 ++++
utilities/recon_cartesian_2d.cpp | 150 ++++++++++++++++++++++++++++++++++++
utilities/test_recon_dataset.cpp | 159 ---------------------------------------
3 files changed, 166 insertions(+), 159 deletions(-)
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index c5e659d..e1b20ff 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -29,6 +29,22 @@ if(FFTW3_FOUND)
install(TARGETS ismrmrd_generate_cartesian_shepp_logan
DESTINATION bin)
+ # Shepp-Logan phantom
+ add_executable(ismrmrd_recon_cartesian_2d
+ recon_cartesian_2d.cpp)
+ if(WIN32)
+ target_link_libraries( ismrmrd_recon_cartesian_2d
+ ismrmrd
+ ${FFTW3_LIBRARIES})
+ else(WIN32)
+ target_link_libraries( ismrmrd_recon_cartesian_2d
+ ismrmrd
+ ${Boost_PROGRAM_OPTIONS_LIBRARY}
+ ${FFTW3_LIBRARIES})
+ endif(WIN32)
+ install(TARGETS ismrmrd_recon_cartesian_2d
+ DESTINATION bin)
+
else(FFTW3_FOUND)
message("FFTW3 NOT Found, cannot build utilities")
endif(FFTW3_FOUND)
diff --git a/utilities/recon_cartesian_2d.cpp b/utilities/recon_cartesian_2d.cpp
new file mode 100644
index 0000000..f6b45d0
--- /dev/null
+++ b/utilities/recon_cartesian_2d.cpp
@@ -0,0 +1,150 @@
+/*
+ * test_recon_dataset.cpp
+ *
+ * Created on: Sep 6, 2012
+ * Author: Michael S. Hansen (michael.hansen at nih.gov)
+ *
+ */
+
+#include <iostream>
+#include "ismrmrd/ismrmrd.h"
+#include "ismrmrd/dataset.h"
+#include "ismrmrd/xml.h"
+#include "fftw3.h"
+
+//Helper function for the FFTW library
+void circshift(complex_float_t *out, const complex_float_t *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))
+
+void print_usage(const char* application)
+{
+ std::cout << "Usage:" << std::endl;
+ std::cout << " - " << application << " <HDF5_FILENAME> " << std::endl;
+}
+
+/* MAIN APPLICATION */
+int main(int argc, char** argv)
+{
+ if (argc < 2) {
+ print_usage(argv[0]);
+ return -1;
+ }
+
+ std::string datafile(argv[1]);
+
+ std::cout << "Simple ISMRMRD Reconstruction program" << std::endl;
+ std::cout << " - filename: " << datafile << std::endl;
+
+ //Let's open the existing dataset
+ ISMRMRD::Dataset d(datafile.c_str(),"dataset", false);
+
+ std::string xml;
+ d.readHeader(xml);
+ ISMRMRD::IsmrmrdHeader hdr;
+ ISMRMRD::deserialize(xml.c_str(),hdr);
+
+ //Let's print some information from the header
+ if (hdr.encoding.size() != 1) {
+ std::cout << "Number of encoding spaces: " << hdr.encoding.size() << std::endl;
+ std::cout << "This simple reconstruction application only supports one encoding space" << std::endl;
+ return -1;
+ }
+
+ ISMRMRD::EncodingSpace e_space = hdr.encoding[0].encodedSpace;
+ ISMRMRD::EncodingSpace r_space = hdr.encoding[0].reconSpace;
+
+ std::cout << "Encoding Matrix Size : [" << e_space.matrixSize.x << ", " << e_space.matrixSize.y << ", " << e_space.matrixSize.z << "]" << std::endl;
+ std::cout << "Reconstruction Matrix Size : [" << r_space.matrixSize.x << ", " << r_space.matrixSize.y << ", " << r_space.matrixSize.z << "]" << std::endl;
+ std::cout << "Number of acquisitions : " << d.getNumberOfAcquisitions() << std::endl;
+
+ if (e_space.matrixSize.z != 1) {
+ std::cout << "This simple reconstruction application only supports 2D encoding spaces" << std::endl;
+ return -1;
+ }
+
+ //Allocate a buffer for the data
+ std::vector<size_t> dims;
+ dims.push_back(e_space.matrixSize.x);
+ dims.push_back(e_space.matrixSize.y);
+ ISMRMRD::NDArray<complex_float_t> buffer(dims);
+
+ //Now loop through and copy data
+ unsigned int number_of_acquisitions = d.getNumberOfAcquisitions();
+ ISMRMRD::Acquisition acq;
+ for (unsigned int i = 0; i < number_of_acquisitions; i++) {
+ //Read one acquisition at a time
+ d.readAcquisition(i, acq);
+
+ //Copy data, we should probably be more careful here and do more tests....
+ //We are not considering multiple channels here.
+ unsigned int offset = acq.idx().kspace_encode_step_1*dims[0];
+ memcpy(&buffer.getData()[offset], acq.getData(),sizeof(complex_float_t)*dims[0]);
+ }
+
+ //Let's FFT the k-space to image
+ fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*buffer.getNumberOfElements());
+
+ if (!tmp) {
+ std::cout << "Error allocating temporary storage for FFTW" << std::endl;
+ return -1;
+ }
+
+ //FFTSHIFT
+ fftshift(reinterpret_cast<complex_float_t*>(tmp), buffer.getData(), dims[0], dims[1]);
+
+ //Create the FFTW plan
+ fftwf_plan p = fftwf_plan_dft_2d(dims[1], dims[0], tmp ,tmp, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+ //Execute the FFT
+ fftwf_execute(p);
+
+ //FFTSHIFT
+ fftshift( buffer.getData(), reinterpret_cast<std::complex<float>*>(tmp), dims[0], dims[1]);
+
+ //Clean up.
+ fftwf_destroy_plan(p);
+ fftwf_free(tmp);
+
+ //Allocate an image
+ ISMRMRD::Image img_out;
+ img_out.data_type(ISMRMRD::ISMRMRD_FLOAT);
+ uint16_t matrix_size[3];
+ matrix_size[0] = r_space.matrixSize.x;
+ matrix_size[1] = r_space.matrixSize.y;
+ matrix_size[2] = 1;
+ img_out.matrix_size(matrix_size);
+ img_out.channels(1);
+
+ //f there is oversampling in the readout direction remove it
+ //Take the magnitude
+ size_t offset = ((e_space.matrixSize.x - r_space.matrixSize.x)>>1);
+ for (unsigned int y = 0; y < matrix_size[1]; y++) {
+ for (unsigned int x = 0; x < matrix_size[0]; x++) {
+ static_cast<float*>(img_out.getData())[y*matrix_size[0]+x] =
+ std::abs(buffer.getData()[y*dims[0] + x + offset]);
+ }
+ }
+
+ // The following are extra guidance we can put in the image header
+ img_out.image_type() = ISMRMRD::ISMRMRD_IMTYPE_MAGNITUDE;
+ img_out.slice() = 0;
+ img_out.field_of_view()[0] = r_space.fieldOfView_mm.x;
+ img_out.field_of_view()[1] = r_space.fieldOfView_mm.y;
+ img_out.field_of_view()[2] = r_space.fieldOfView_mm.z;
+ //And so on
+
+ //Let's write the reconstructed image into the same data file
+ d.appendImage("myimage", ISMRMRD::ISMRMRD_BLOCKMODE_ARRAY, img_out);
+
+ return 0;
+}
diff --git a/utilities/test_recon_dataset.cpp b/utilities/test_recon_dataset.cpp
deleted file mode 100644
index 390e70e..0000000
--- a/utilities/test_recon_dataset.cpp
+++ /dev/null
@@ -1,159 +0,0 @@
-/*
- * test_recon_dataset.cpp
- *
- * Created on: Sep 6, 2012
- * Author: Michael S. Hansen (michael.hansen at nih.gov)
- *
- */
-
-#include <iostream>
-#include "ismrmrd.h"
-#include "ismrmrd.hxx"
-#include "ismrmrd_hdf5.h"
-#include "fftw3.h"
-
-//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))
-
-void print_usage(const char* application)
-{
- std::cout << "Usage:" << std::endl;
- std::cout << " - " << application << " <HDF5_FILENAME> " << "<SCHEMA_FILENAME>" << std::endl;
-}
-
-/* MAIN APPLICATION */
-int main(int argc, char** argv)
-{
- if (argc < 3) {
- print_usage(argv[0]);
- return -1;
- }
-
- std::string datafile(argv[1]);
- std::string schemafile(argv[2]);
-
- std::cout << "Simple ISMRMRD Reconstruction program" << std::endl;
- std::cout << " - filename: " << datafile << std::endl;
- std::cout << " - schema : " << schemafile << std::endl;
-
- //Let's open the dataset
- ISMRMRD::IsmrmrdDataset d(datafile.c_str(),"dataset");
-
- //We will start by getting the header and turning it into a C++ class
- //In order to do true validation of the XML, we will use the XML schema
- xml_schema::properties props;
- props.schema_location ("http://www.ismrm.org/ISMRMRD", schemafile);
-
- boost::shared_ptr<std::string> xml = d.readHeader();
- std::istringstream str_stream(*xml, std::stringstream::in);
-
- boost::shared_ptr<ISMRMRD::ismrmrdHeader> cfg;
-
- try {
- cfg = boost::shared_ptr<ISMRMRD::ismrmrdHeader>(ISMRMRD::ismrmrdHeader_ (str_stream,0,props));
- } catch (const xml_schema::exception& e) {
- std::cout << "Failed to parse XML Parameters: " << e.what() << std::endl;
- }
-
- //Let's print some information from the header
- ISMRMRD::ismrmrdHeader::encoding_sequence e_seq = cfg->encoding();
- if (e_seq.size() != 1) {
- std::cout << "Number of encoding spaces: " << e_seq.size() << std::endl;
- std::cout << "This simple reconstruction application only supports one encoding space" << std::endl;
- return -1;
- }
-
- ISMRMRD::encodingSpaceType e_space = (*e_seq.begin()).encodedSpace();
- ISMRMRD::encodingSpaceType r_space = (*e_seq.begin()).reconSpace();
- ISMRMRD::encodingLimitsType e_limits = (*e_seq.begin()).encodingLimits();
-
- std::cout << "Encoding Matrix Size : [" << e_space.matrixSize().x() << ", " << e_space.matrixSize().y() << ", " << e_space.matrixSize().z() << "]" << std::endl;
- std::cout << "Reconstruction Matrix Size : [" << r_space.matrixSize().x() << ", " << r_space.matrixSize().y() << ", " << r_space.matrixSize().z() << "]" << std::endl;
- std::cout << "Number of acquisitions : " << d.getNumberOfAcquisitions() << std::endl;
-
- if (e_space.matrixSize().z() != 1) {
- std::cout << "This simple reconstruction application only supports 2D encoding spaces" << std::endl;
- return -1;
- }
-
-
- //Allocate a buffer for the data
- ISMRMRD::NDArrayContainer< std::complex<float> > buffer;
- buffer.dimensions_.push_back(e_space.matrixSize().x());
- buffer.dimensions_.push_back(e_space.matrixSize().y());
- buffer.resize(e_space.matrixSize().x()*e_space.matrixSize().y(), std::complex<float>(0.0,0.0));
-
-
- //Now loop through and copy data
- unsigned int number_of_acquisitions = d.getNumberOfAcquisitions();
- for (unsigned int i = 0; i < number_of_acquisitions; i++) {
- //Read one acquisition at a time
- boost::shared_ptr<ISMRMRD::Acquisition> acq = d.readAcquisition(i);
-
- //Copy data, we should probably be more careful here and do more tests....
- //We are not considering multiple channels here.
- unsigned int offset = acq->getIdx().kspace_encode_step_1*buffer.dimensions_[0];
- memcpy(&buffer[offset],&acq->getData()[0],sizeof(float)*2*buffer.dimensions_[0]);
- }
-
- //Let's FFT the k-space to image
- fftwf_complex* tmp = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*buffer.data_.size());
-
- if (!tmp) {
- std::cout << "Error allocating temporary storage for FFTW" << std::endl;
- return -1;
- }
-
- fftshift(reinterpret_cast<std::complex<float>*>(tmp),&buffer.data_[0],buffer.dimensions_[0],buffer.dimensions_[1]);
-
- //Create the FFTW plan
- fftwf_plan p = fftwf_plan_dft_2d(buffer.dimensions_[1], buffer.dimensions_[0], tmp,tmp, FFTW_BACKWARD, FFTW_ESTIMATE);
-
- fftwf_execute(p);
-
- fftshift(&buffer.data_[0],reinterpret_cast<std::complex<float>*>(tmp),buffer.dimensions_[0],buffer.dimensions_[1]);
-
- //Clean up.
- fftwf_destroy_plan(p);
- fftwf_free(tmp);
-
- //Now, let's remove the oversampling in the readout and take the magnitude
- //Allocate a buffer for the data
- ISMRMRD::NDArrayContainer< float > img;
- img.dimensions_.push_back(r_space.matrixSize().x());
- img.dimensions_.push_back(r_space.matrixSize().y());
- img.data_.resize(r_space.matrixSize().x()*r_space.matrixSize().y(), 0.0);
-
- for (unsigned int y = 0; y < img.dimensions_[1]; y++) {
- for (unsigned int x = 0; x < img.dimensions_[0]; x++) {
- img.data_[y*img.dimensions_[0]+x] =
- std::abs(buffer.data_[y*buffer.dimensions_[0] + x + ((e_space.matrixSize().x()-r_space.matrixSize().x())>>1)]);
- }
- }
-
- //Let's write the reconstructed image straight in the same data file
- ISMRMRD::ImageHeader img_h;
- img_h.channels = 1;
- img_h.image_data_type = ISMRMRD::DATA_COMPLEX_FLOAT; //This is actually just guidance
- img_h.image_type = ISMRMRD::TYPE_COMPLEX; //This is actually just guidance
- img_h.slice = 0;
- //And so on
-
- //Now append, we will append image and header separately (in two different datasets)
- d.appendImageHeader(img_h,"myimage.head");
- d.appendArray(img,"myimage.img");
-
-
- 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