[ismrmrd] 09/281: Added examples 2D/3D and spiral with Matlab code
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:00:49 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 3fbc194376edf877c652a7d09772dc7462b7868a
Author: Michael S. Hansen <michael.hansen at nih.gov>
Date: Thu Aug 9 20:38:52 2012 -0400
Added examples 2D/3D and spiral with Matlab code
---
examples/matlab/ismrmrd_header2struct.m | 25 +++++
examples/matlab/simple_cartesian_recon.m | 82 ++++++++++++++
examples/matlab/simple_gridder.m | 100 +++++++++++++++++
examples/matlab/simple_spiral_recon.m | 78 +++++++++++++
examples/matlab/xml2struct.m | 183 +++++++++++++++++++++++++++++++
schema/ismrmrd.xsd | 11 +-
6 files changed, 478 insertions(+), 1 deletion(-)
diff --git a/examples/matlab/ismrmrd_header2struct.m b/examples/matlab/ismrmrd_header2struct.m
new file mode 100644
index 0000000..34e2370
--- /dev/null
+++ b/examples/matlab/ismrmrd_header2struct.m
@@ -0,0 +1,25 @@
+function s = ismrmrd_header2struct(header)
+%
+%
+% s = ismrmrd_header2struct(header)
+%
+% Converts and ISMRMRD xml header to a struct.
+%
+
+ if (iscell(header)),
+ header_string = header{1};
+ else
+ header_string = header;
+ end
+
+ if (ischar(header_string) < 1),
+ error('Malformed input header. Is not a character string.');
+ end
+
+ tmp_nam = tempname;
+ f=fopen(tmp_nam,'w');
+ fwrite(f,header_string);
+ fclose(f);
+
+ s = xml2struct(tmp_nam);
+return
\ No newline at end of file
diff --git a/examples/matlab/simple_cartesian_recon.m b/examples/matlab/simple_cartesian_recon.m
new file mode 100644
index 0000000..41fe0ba
--- /dev/null
+++ b/examples/matlab/simple_cartesian_recon.m
@@ -0,0 +1,82 @@
+function [images] = simple_cartesian_recon(ismrmrdfile)
+%
+% [images] = simple_cartesian_recon(ismrmrdfile)
+%
+% Simple example of how to reconstruct a 2D/3D cartesian dataset stored
+% in ISMRMRD dataformat. Reconstruction handles partial fourier,
+% oversampling, etc.
+%
+% Michael S. Hansen (michael.hansen at nih.gov), 2012
+%
+
+
+%read the header
+header = ismrmrd_header2struct(h5read(ismrmrdfile,'/dataset/xml'));
+
+
+%Is this a spiral acquisition
+if (~strcmp(header.ismrmrdHeader.encoding.trajectory.Text,'cartesian')),
+ error('This is not a cartesian dataset');
+end
+
+%if (str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.z.Text) <= 1),
+% error('This is not a 3D dataset');
+%end
+
+%Let's get the matrix size
+matrix_size = [str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.x.Text), ...
+ str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.y.Text), ...
+ str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.z.Text)];
+
+
+recon_size = [str2num(header.ismrmrdHeader.encoding.reconSpace.matrixSize.x.Text), ...
+ str2num(header.ismrmrdHeader.encoding.reconSpace.matrixSize.y.Text), ...
+ str2num(header.ismrmrdHeader.encoding.reconSpace.matrixSize.z.Text)];
+
+
+%Let's load the data
+raw_data = h5read(ismrmrdfile,'/dataset/data');
+
+channels = max(raw_data.head.active_channels);
+samples = max(raw_data.head.number_of_samples);
+center_line = str2num(header.ismrmrdHeader.encoding.encodingLimits.kspace_encoding_step_1.center.Text);
+center_partition = str2num(header.ismrmrdHeader.encoding.encodingLimits.kspace_encoding_step_2.center.Text);
+
+if (samples ~= matrix_size(1)),
+ error('Mismatch between number of samples and matrix size');
+end
+
+buffer = zeros(matrix_size(1),matrix_size(2),matrix_size(3),channels);
+counter = 0;
+for p=1:length(raw_data.head.flags),
+ if (bitget(uint64(raw_data.head.flags(p)),19)), %if this is noise, we will skip it
+ continue;
+ end
+ line_number = (int32(raw_data.head.idx.kspace_encode_step_1(p))-center_line)+bitshift(matrix_size(2),-1);
+ partition_number = (int32(raw_data.head.idx.kspace_encode_step_2(p))-center_partition)+bitshift(matrix_size(3),-1);
+ buffer(:,line_number+1,partition_number+1,:) = reshape(complex(raw_data.data{p}.real, raw_data.data{p}.imag), samples, 1, 1, channels);
+
+ if (bitget(uint64(raw_data.head.flags(p)),8)), %Is this the last in slice? We should make an image
+ fprintf('Reconstructing image %d....', counter+1);
+
+ img = fftshift(ifft(ifftshift(buffer),[],1)) ./ sqrt(size(buffer,1));
+ img = fftshift(ifft(ifftshift(img),[],2)) ./ sqrt(size(img,2));
+ if (matrix_size(3) > 1),
+ img = fftshift(ifft(ifftshift(img),[],3)) ./ sqrt(size(img,3));
+ end
+
+ img = img((1:recon_size(1))+bitshift(matrix_size(1)-recon_size(1),-1), ...
+ (1:recon_size(2))+bitshift(matrix_size(2)-recon_size(2),-1), ...
+ (1:recon_size(3))+bitshift(matrix_size(3)-recon_size(3),-1), :);
+
+ images(:,:,:,counter+1) = sqrt(sum(abs(img).^2,4)); %RMS coil combination
+ counter = counter + 1;
+ fprintf('done\n');
+ buffer = zeros(matrix_size(1),matrix_size(2),matrix_size(3), channels);
+ end
+end
+
+%Let's just show the middle slice image
+imagesc(images(:,:,bitshift(size(images,3),-1)+1, 1));colormap(gray);axis image;
+
+
diff --git a/examples/matlab/simple_gridder.m b/examples/matlab/simple_gridder.m
new file mode 100644
index 0000000..cb65fb5
--- /dev/null
+++ b/examples/matlab/simple_gridder.m
@@ -0,0 +1,100 @@
+function img = simple_gridder(co,data,weight,matrix_size)
+%
+% img = simple_gridder(co,data,weight,matrix_size)
+%
+% A very simple gridder written purely in Matlab.
+% Only tested for 2D
+%
+% Input arguments:
+% - co [N,2], k-space coordinates in the range [-0.5:0.5]
+% - data, vector of N complex data points
+% - weight, vector of N data weights (density compensation)
+% - matrix_size, 2-element vector, e.g. [128 128]
+%
+% Michael S. Hansen (michael.hansen at nih.gov), 2012
+%
+%
+
+over_sampling = 2.0;
+kernel_width = 8;
+kernel_samples = 100000;
+kernel_beta = 18.5547;
+matrix_size_oversampled = matrix_size*over_sampling;
+
+kernel_table = 1.0-(2.0*[-bitshift(kernel_samples,-1):(bitshift(kernel_samples,-1)-1)]/kernel_samples).^2;
+
+kernel_table(kernel_table<0) = 0;
+kernel_table=sqrt(kernel_table);
+kernel_table = bessi0(kernel_beta .* kernel_table);
+knorm = sum(kernel_table(:))/numel(kernel_table);
+kernel_table = kernel_table ./ knorm;
+
+grid_cells = ceil(co .* repmat(matrix_size_oversampled,size(co,1),1)-(kernel_width/2));
+kernel_idx = grid_cells-(co .* repmat(matrix_size_oversampled,size(co,1),1));
+grid_cells = uint32(grid_cells+repmat(bitshift(matrix_size_oversampled+kernel_width,-1),size(co,1),1));
+kernel_idx = uint32(floor(((kernel_idx + kernel_width/2)/kernel_width)*kernel_samples));
+kernel_step = uint32(floor(kernel_samples/kernel_width));
+
+%Calculate deapodization
+kern = -bitshift(kernel_width,-1):(-bitshift(kernel_width,-1)+kernel_width-1);
+kern = 1.0-(2.0*kern/kernel_width).^2;
+kern = bessi0(kernel_beta .* kern);
+kern = repmat(kern,kernel_width,1) .* repmat(kern',1,kernel_width);
+deapodization = zeros(matrix_size_oversampled);
+deapodization((1:kernel_width) + bitshift(matrix_size_oversampled(1)-kernel_width,-1), ...
+ (1:kernel_width) + bitshift(matrix_size_oversampled(2)-kernel_width,-1)) = kern;
+deapodization = fftshift(ifftn(ifftshift(deapodization))) ./ sqrt(numel(deapodization));
+
+%Do convolution
+oversampled_grid_padded = zeros(matrix_size_oversampled+kernel_width);
+for y=1:(kernel_width),
+ for x=1:(kernel_width),
+ oversampled_grid_padded = oversampled_grid_padded + ...
+ accumarray([grid_cells(:,1)+(x-1),grid_cells(:,2)+(y-1)], ...
+ weight(:).*data(:).*(kernel_table(kernel_idx(:,1)+(x-1)*kernel_step+1).*kernel_table(kernel_idx(:,2)+(y-1).*kernel_step+1))', ...
+ matrix_size_oversampled+kernel_width);
+ end
+end
+
+osps = size(oversampled_grid_padded);
+oss = matrix_size_oversampled;
+
+%Let's just get rid of the padding, it should really be folded in
+oversampled_grid = oversampled_grid_padded(bitshift(osps(1)-oss(1),-1):bitshift(osps(1)-oss(1),-1)+matrix_size_oversampled(1)-1, ...
+ bitshift(osps(2)-oss(2),-1):bitshift(osps(1)-oss(2),-1)+matrix_size_oversampled(2)-1);
+
+%FFT to image space
+img = fftshift(ifftn(ifftshift(oversampled_grid))) ./ sqrt(numel(oversampled_grid));
+
+%Deapodization
+img = img ./ deapodization;
+
+%Remove oversampling
+img = img((1:matrix_size(1))+bitshift(matrix_size_oversampled(1)-matrix_size(1),-1), ...
+ (1:matrix_size(2))+bitshift(matrix_size_oversampled(2)-matrix_size(2),-1));
+
+%TODO
+% - fold padding back in
+
+return
+
+function ans = bessi0(x)
+
+ax = abs(x);
+ans = zeros(size(x));
+
+% ax<3.75
+k = find(ax<3.75);
+y=x(k)./3.75;
+y=y.^2;
+ans(k)=1.0+y.*(3.5156229+y.*(3.0899424+y.*(1.2067492 ...
+ +y.*(0.2659732+y.*(0.360768e-1+y.*0.45813e-2)))));
+
+% ax>=3.75
+k = find(ax>=3.75);
+y=3.75./ax(k);
+ans(k)=(exp(ax(k))./sqrt(ax(k))).*(0.39894228+y.*(0.1328592e-1 ...
+ +y.*(0.225319e-2+y.*(-0.157565e-2+y.*(0.916281e-2 ...
+ +y.*(-0.2057706e-1+y.*(0.2635537e-1+y.*(-0.1647633e-1 ...
+ +y.*0.392377e-2))))))));
+return
\ No newline at end of file
diff --git a/examples/matlab/simple_spiral_recon.m b/examples/matlab/simple_spiral_recon.m
new file mode 100644
index 0000000..0851211
--- /dev/null
+++ b/examples/matlab/simple_spiral_recon.m
@@ -0,0 +1,78 @@
+function [images] = simple_spiral_recon(ismrmrdfile)
+%
+% [images] = simple_spiral_recon(ismrmrdfile)
+%
+% Simple example of how to reconstruct a 2D spiral dataset stored
+% in ISMRMRD dataformat
+%
+% Michael S. Hansen (michael.hansen at nih.gov), 2012
+%
+
+header = ismrmrd_header2struct(h5read(ismrmrdfile,'/dataset/xml'));
+
+%Is this a spiral acquisition
+if (~strcmp(header.ismrmrdHeader.encoding.trajectory.Text,'spiral')),
+ error('This is not a spiral dataset');
+end
+
+%Let's get the matrix size
+matrix_size = [str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.x.Text), ...
+ str2num(header.ismrmrdHeader.encoding.encodedSpace.matrixSize.y.Text)];
+
+%Let's load the data
+raw_data = h5read(ismrmrdfile,'/dataset/data');
+
+interleaves = max(raw_data.head.idx.kspace_encode_step_1)+1;
+repetitions = max(raw_data.head.idx.repetition)+1;
+samples = max(raw_data.head.number_of_samples);
+samples_to_skip_end = max(raw_data.head.discard_post);
+samples_to_skip_start = max(raw_data.head.discard_pre);
+net_samples = samples-samples_to_skip_end-samples_to_skip_start;
+channels = max(raw_data.head.active_channels);
+
+trajectory = zeros(2,net_samples,interleaves);
+data = zeros(net_samples,interleaves,channels);
+
+images = [];
+
+counter = 0;
+for p=1:length(raw_data.head.flags),
+
+ if (bitget(uint64(raw_data.head.flags(p)),19)), %if this is noise, we will skip it
+ continue;
+ end
+
+ d = reshape(complex(raw_data.data{p}.real, raw_data.data{p}.imag), samples, channels);
+ t = reshape(raw_data.traj{p}, raw_data.head.trajectory_dimensions(p), samples);
+ current_interleave = raw_data.head.idx.kspace_encode_step_1(p)+1;
+ start_sample = samples_to_skip_start+1;
+ end_sample = samples-samples_to_skip_end;
+
+ data(:,current_interleave,:) = reshape(d(start_sample:end_sample,:), net_samples, 1, channels);
+ trajectory(:,:,current_interleave) = t(:,start_sample:end_sample);
+
+ if (bitget(uint64(raw_data.head.flags(p)),8)), %Is this the last in slice? We should make an image
+ fprintf('Reconstructing image %d....', counter+1);
+ co = permute(trajectory(:,:),[2 1]);
+
+ %Some simple density compensation
+ k = complex(co(:,1),co(:,2));
+ g = complex(diff(co(:,1)),diff(co(:,2)));
+ g(end+1) = g(end);
+ weights = abs(g(:)) .* abs(sin(angle(g(:))-angle(k(:)))); %Estimating weights from Meyer et al. Magn Reson Med. 1992 Dec;28(2):202-13.
+
+ %Loop over coils and grid images
+ for c=1:size(data,3),
+ img(:,:,c) = simple_gridder(co,data(:,:,c),weights, matrix_size);
+ end;
+ images(:,:,counter+1) = sqrt(sum(abs(img).^2,3)); %RMS coil combination
+ counter = counter + 1;
+ fprintf('done\n');
+ end
+
+end
+
+%Let's just display the first image
+imagesc(images(:,:,1));colormap(gray);axis image;
+
+return
diff --git a/examples/matlab/xml2struct.m b/examples/matlab/xml2struct.m
new file mode 100644
index 0000000..dcd85a2
--- /dev/null
+++ b/examples/matlab/xml2struct.m
@@ -0,0 +1,183 @@
+function [ s ] = xml2struct( file )
+%Convert xml file into a MATLAB structure
+% [ s ] = xml2struct( file )
+%
+% A file containing:
+% <XMLname attrib1="Some value">
+% <Element>Some text</Element>
+% <DifferentElement attrib2="2">Some more text</Element>
+% <DifferentElement attrib3="2" attrib4="1">Even more text</DifferentElement>
+% </XMLname>
+%
+% Will produce:
+% s.XMLname.Attributes.attrib1 = "Some value";
+% s.XMLname.Element.Text = "Some text";
+% s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2";
+% s.XMLname.DifferentElement{1}.Text = "Some more text";
+% s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2";
+% s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1";
+% s.XMLname.DifferentElement{2}.Text = "Even more text";
+%
+% Please note that the following characters are substituted
+% '-' by '_dash_', ':' by '_colon_' and '.' by '_dot_'
+%
+% Written by W. Falkena, ASTI, TUDelft, 21-08-2010
+% Attribute parsing speed increased by 40% by A. Wanner, 14-6-2011
+% Added CDATA support by I. Smirnov, 20-3-2012
+%
+% Modified by X. Mo, University of Wisconsin, 12-5-2012
+
+ if (nargin < 1)
+ clc;
+ help xml2struct
+ return
+ end
+
+ if isa(file, 'org.apache.xerces.dom.DeferredDocumentImpl') || isa(file, 'org.apache.xerces.dom.DeferredElementImpl')
+ % input is a java xml object
+ xDoc = file;
+ else
+ %check for existance
+ if (exist(file,'file') == 0)
+ %Perhaps the xml extension was omitted from the file name. Add the
+ %extension and try again.
+ if (isempty(strfind(file,'.xml')))
+ file = [file '.xml'];
+ end
+
+ if (exist(file,'file') == 0)
+ error(['The file ' file ' could not be found']);
+ end
+ end
+ %read the xml file
+ xDoc = xmlread(file);
+ end
+
+ %parse xDoc into a MATLAB structure
+ s = parseChildNodes(xDoc);
+
+end
+
+% ----- Subfunction parseChildNodes -----
+function [children,ptext,textflag] = parseChildNodes(theNode)
+ % Recurse over node children.
+ children = struct;
+ ptext = struct; textflag = 'Text';
+ if hasChildNodes(theNode)
+ childNodes = getChildNodes(theNode);
+ numChildNodes = getLength(childNodes);
+
+ for count = 1:numChildNodes
+ theChild = item(childNodes,count-1);
+ [text,name,attr,childs,textflag] = getNodeData(theChild);
+
+ if (~strcmp(name,'#text') && ~strcmp(name,'#comment') && ~strcmp(name,'#cdata_dash_section'))
+ %XML allows the same elements to be defined multiple times,
+ %put each in a different cell
+ if (isfield(children,name))
+ if (~iscell(children.(name)))
+ %put existsing element into cell format
+ children.(name) = {children.(name)};
+ end
+ index = length(children.(name))+1;
+ %add new element
+ children.(name){index} = childs;
+ if(~isempty(fieldnames(text)))
+ children.(name){index} = text;
+ end
+ if(~isempty(attr))
+ children.(name){index}.('Attributes') = attr;
+ end
+ else
+ %add previously unknown (new) element to the structure
+ children.(name) = childs;
+ if(~isempty(text) && ~isempty(fieldnames(text)))
+ children.(name) = text;
+ end
+ if(~isempty(attr))
+ children.(name).('Attributes') = attr;
+ end
+ end
+ else
+ ptextflag = 'Text';
+ if (strcmp(name, '#cdata_dash_section'))
+ ptextflag = 'CDATA';
+ elseif (strcmp(name, '#comment'))
+ ptextflag = 'Comment';
+ end
+
+ %this is the text in an element (i.e., the parentNode)
+ if (~isempty(regexprep(text.(textflag),'[\s]*','')))
+ if (~isfield(ptext,ptextflag) || isempty(ptext.(ptextflag)))
+ ptext.(ptextflag) = text.(textflag);
+ else
+ %what to do when element data is as follows:
+ %<element>Text <!--Comment--> More text</element>
+
+ %put the text in different cells:
+ % if (~iscell(ptext)) ptext = {ptext}; end
+ % ptext{length(ptext)+1} = text;
+
+ %just append the text
+ ptext.(ptextflag) = [ptext.(ptextflag) text.(textflag)];
+ end
+ end
+ end
+
+ end
+ end
+end
+
+% ----- Subfunction getNodeData -----
+function [text,name,attr,childs,textflag] = getNodeData(theNode)
+ % Create structure of node info.
+
+ %make sure name is allowed as structure name
+ name = toCharArray(getNodeName(theNode))';
+ name = strrep(name, '-', '_dash_');
+ name = strrep(name, ':', '_colon_');
+ name = strrep(name, '.', '_dot_');
+
+ attr = parseAttributes(theNode);
+ if (isempty(fieldnames(attr)))
+ attr = [];
+ end
+
+ %parse child nodes
+ [childs,text,textflag] = parseChildNodes(theNode);
+
+ if (isempty(fieldnames(childs)) && isempty(fieldnames(text)))
+ %get the data of any childless nodes
+ % faster than if any(strcmp(methods(theNode), 'getData'))
+ % no need to try-catch (?)
+ % faster than text = char(getData(theNode));
+ text.(textflag) = toCharArray(getTextContent(theNode))';
+ end
+
+end
+
+% ----- Subfunction parseAttributes -----
+function attributes = parseAttributes(theNode)
+ % Create attributes structure.
+
+ attributes = struct;
+ if hasAttributes(theNode)
+ theAttributes = getAttributes(theNode);
+ numAttributes = getLength(theAttributes);
+
+ for count = 1:numAttributes
+ %attrib = item(theAttributes,count-1);
+ %attr_name = regexprep(char(getName(attrib)),'[-:.]','_');
+ %attributes.(attr_name) = char(getValue(attrib));
+
+ %Suggestion of Adrian Wanner
+ str = toCharArray(toString(item(theAttributes,count-1)))';
+ k = strfind(str,'=');
+ attr_name = str(1:(k(1)-1));
+ attr_name = strrep(attr_name, '-', '_dash_');
+ attr_name = strrep(attr_name, ':', '_colon_');
+ attr_name = strrep(attr_name, '.', '_dot_');
+ attributes.(attr_name) = str((k(1)+2):(end-1));
+ end
+ end
+end
\ No newline at end of file
diff --git a/schema/ismrmrd.xsd b/schema/ismrmrd.xsd
index cc9fbf6..4633b55 100644
--- a/schema/ismrmrd.xsd
+++ b/schema/ismrmrd.xsd
@@ -18,6 +18,7 @@
</xs:all>
</xs:complexType>
</xs:element>
+ <xs:element maxOccurs="1" minOccurs="0" name="sequenceTiming" type="sequenceTimingType"/>
<xs:element maxOccurs="1" minOccurs="0" name="userParameters">
<xs:complexType>
<xs:sequence>
@@ -48,7 +49,7 @@
<xs:element minOccurs="0" name="systemVendor" type="xs:string"/>
<xs:element minOccurs="0" name="systemModel" type="xs:string"/>
<xs:element minOccurs="0" name="systemFieldStrength_T" type="xs:float"/>
- <xs:element minOccurs="0" name="receivedChannels" type="xs:unsignedShort"/>
+ <xs:element minOccurs="0" name="receiverChannels" type="xs:unsignedShort"/>
</xs:all>
</xs:complexType>
@@ -118,6 +119,14 @@
</xs:sequence>
</xs:complexType>
+ <xs:complexType name="sequenceTimingType">
+ <xs:sequence>
+ <xs:element minOccurs="1" maxOccurs="unbounded" type="xs:float" name="TR"/>
+ <xs:element minOccurs="1" maxOccurs="unbounded" type="xs:float" name="TE"/>
+ <xs:element minOccurs="0" maxOccurs="unbounded" type="xs:float" name="TI"/>
+ </xs:sequence>
+ </xs:complexType>
+
<xs:complexType name="userParameterLongType">
<xs:all>
<xs:element name="name" type="xs:string"/>
--
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