[ismrmrd] 190/281: A simple matlab reconstruction example.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:01:13 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 a203d7d312403bb8776625bdc6cb1c915640ba30
Author: Souheil Inati <souheil.inati at nih.gov>
Date: Tue Dec 3 13:40:42 2013 -0500
A simple matlab reconstruction example.
---
examples/matlab/example.m | 41 -----------
examples/matlab/recon_simple.m | 158 +++++++++++++++++++++++++++++++++++++++++
2 files changed, 158 insertions(+), 41 deletions(-)
diff --git a/examples/matlab/example.m b/examples/matlab/example.m
deleted file mode 100644
index 84c5051..0000000
--- a/examples/matlab/example.m
+++ /dev/null
@@ -1,41 +0,0 @@
-function example()
-
-filename = 'test.h5';
-
-if (exist(filename, 'file'))
- delete(filename) ;
-end
-
-% Create the dataset
-i = ismrmrd.IsmrmrdDataset(filename, 'dataset');
-
-ah = ismrmrd.AcquisitionHeader();
-ah.center_sample = 42;
-ah.active_channels = 32;
-ah.version = 2;
-ah.position = [3.14, 3.14, 3.14]
-
-a = ismrmrd.Acquisition();
-a.head = ah;
-a.data = (1:512);
-a.traj = (1:3);
-
-i.appendAcquisition(a);
-
-% Write the dataset
-i.close();
-
-% Open the dataset
-i = ismrmrd.IsmrmrdDataset(filename, 'dataset')
-
-% Read the dataset
-nacq = i.getNumberOfAcquisitions();
-a = i.readAcquisition(1);
-
-h = a.head;
-
-h.position
-
-i.close();
-
-end
diff --git a/examples/matlab/recon_simple.m b/examples/matlab/recon_simple.m
new file mode 100644
index 0000000..3cb0a16
--- /dev/null
+++ b/examples/matlab/recon_simple.m
@@ -0,0 +1,158 @@
+%% Working with an existing ISMRMRD data set
+
+% This is a simple example of how to reconstruct images from data
+% acquired on a fully sampled cartesian grid
+%
+% Capabilities:
+% 2D/3D
+% multiple slices/slabs
+% multiple contrasts, repetitions
+%
+% Limitations:
+% only works with a single encoded space
+% fully sampled k-space (no partial fourier or undersampling)
+% doesn't handle averages, phases, segments and sets
+% ignores noise scans (no pre-whitening)
+%
+
+% We first create a data set using the example program like this:
+% ismrmrd_generate_cartesian_shepp_logan -r 5 -C -o shepp-logan.h5
+% This will produce the file shepp-logan.h5 containing an ISMRMRD
+% dataset sampled evenly on the k-space grid -128:0.5:127.5 x -128:127
+% (i.e. oversampled by a factor of 2 in the readout direction)
+% with 8 coils, 5 repetitions and a noise level of 0.5
+% with a noise calibration scan at the beginning
+%
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Loading an existing file %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+filename = 'shepp-logan.h5';
+if exist(filename, 'file')
+ dset = ismrmrd.IsmrmrdDataset(filename, 'dataset');
+else
+ error(['File ' filename ' does not exist. Please generate it.'])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Read some fields from the XML header %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% We need to check if optional fields exists before trying to read them
+% Some functions return java.lang.Integer type and we convert to double
+
+%% Encoding and reconstruction information
+% Matrix size
+enc_Nx = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getMatrixSize.getX;
+enc_Ny = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getMatrixSize.getY;
+enc_Nz = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getMatrixSize.getZ;
+rec_Nx = dset.xmlhdr.getEncoding.get(0).getReconSpace.getMatrixSize.getX;
+rec_Ny = dset.xmlhdr.getEncoding.get(0).getReconSpace.getMatrixSize.getY;
+rec_Nz = dset.xmlhdr.getEncoding.get(0).getReconSpace.getMatrixSize.getZ;
+
+% Field of View
+enc_FOVx = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getFieldOfViewMm.getX;
+enc_FOVy = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getFieldOfViewMm.getY;
+enc_FOVz = dset.xmlhdr.getEncoding.get(0).getEncodedSpace.getFieldOfViewMm.getZ;
+rec_FOVx = dset.xmlhdr.getEncoding.get(0).getReconSpace.getFieldOfViewMm.getX;
+rec_FOVy = dset.xmlhdr.getEncoding.get(0).getReconSpace.getFieldOfViewMm.getY;
+rec_FOVz = dset.xmlhdr.getEncoding.get(0).getReconSpace.getFieldOfViewMm.getZ;
+
+% Number of slices
+if isempty(dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getSlice)
+ nSlices = 1;
+else
+ nSlices = dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getSlice.getMaximum;
+end
+
+% Number of coils, in the system information
+nCoils = double(dset.xmlhdr.getAcquisitionSystemInformation.getReceiverChannels);
+
+% Repetitions, contrasts etc.
+if isempty(dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getRepetition)
+ nReps = 1;
+else
+ nReps = double(dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getRepetition.getMaximum);
+end
+
+if isempty(dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getContrast)
+ nContrasts = 1;
+else
+ nContrasts = double(dset.xmlhdr.getEncoding.get(0).getEncodingLimits.getContrast.getMaximum);
+end
+
+% TODO add the other possibilites
+
+%% Read all the data
+% Reading can be done one acquisition (or chunk) at a time,
+% but this is much faster for data sets that fit into RAM.
+D = dset.readAcquisition();
+
+% Note: can select a single acquisition or header from the block, e.g.
+% acq = D.select(5);
+% hdr = D.head.select(5);
+% or you can work with them all together
+
+%% Ignore noise scans
+% TODO add a pre-whitening example
+% Find the first non-noise scan
+% This is how to check if a flag is set in the acquisition header
+isNoise = D.head.flagIsSet(D.head.FLAGS.ACQ_IS_NOISE_MEASUREMENT);
+firstScan = find(isNoise==0,1,'first');
+if firstScan > 1
+ noise = D.select(1:firstScan-1);
+else
+ noise = [];
+end
+meas = D.select(firstScan:D.getNumber);
+clear D;
+
+%% Reconstruct images
+% Since the entire file is in memory we can use random access
+% Loop over repetitions, contrasts, slices
+reconImages = {};
+nimages = 0;
+for rep = 1:nReps
+ for contrast = 1:nContrasts
+ for slice = 1:nSlices
+ % Initialize the K-space storage array
+ K = zeros(enc_Nx, enc_Ny, enc_Nz, nCoils);
+ % Select the appropriate measurements from the data
+ acqs = find( (meas.head.idx.contrast==(contrast-1)) ...
+ & (meas.head.idx.repetition==(rep-1)) ...
+ & (meas.head.idx.slice==(slice-1)));
+ for p = 1:length(acqs)
+ ky = meas.head.idx.kspace_encode_step_1(acqs(p)) + 1;
+ kz = meas.head.idx.kspace_encode_step_2(acqs(p)) + 1;
+ K(:,ky,kz,:) = meas.data{acqs(p)};
+ end
+ % Reconstruct in x
+ K = fftshift(ifft(fftshift(K,1),[],1),1);
+ % Chop if needed
+ if (enc_Nx == rec_Nx)
+ im = K;
+ else
+ ind1 = floor((enc_Nx - rec_Nx)/2)+1;
+ ind2 = floor((enc_Nx - rec_Nx)/2)+rec_Nx;
+ im = K(ind1:ind2,:,:,:);
+ end
+ % Reconstruct in y then z
+ im = fftshift(ifft(fftshift(im,2),[],2),2);
+ if size(im,3)>1
+ im = fftshift(ifft(fftshift(im,3),[],3),3);
+ end
+
+ % Combine SOS across coils
+ im = sqrt(sum(abs(im).^2,4));
+
+ % Append
+ nimages = nimages + 1;
+ reconImages{nimages} = im;
+ end
+ end
+end
+
+%% Display the first image
+figure
+colormap gray
+imagesc(reconImages{1}); axis image; axis off; colorbar;
--
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