[ismrmrd] 164/177: Revived Michael's matlab spiral recon example
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:02:15 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 e8849ff7ac10db3a54b634aede144523c05f4a15
Author: Souheil Inati <souheil.inati at nih.gov>
Date: Fri Oct 24 15:41:22 2014 -0400
Revived Michael's matlab spiral recon example
---
examples/matlab/old/simple_gridder.m | 100 ++++++++++++++++++++++++++++++
examples/matlab/old/simple_spiral_recon.m | 78 +++++++++++++++++++++++
2 files changed, 178 insertions(+)
diff --git a/examples/matlab/old/simple_gridder.m b/examples/matlab/old/simple_gridder.m
new file mode 100644
index 0000000..cb65fb5
--- /dev/null
+++ b/examples/matlab/old/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/old/simple_spiral_recon.m b/examples/matlab/old/simple_spiral_recon.m
new file mode 100644
index 0000000..3815580
--- /dev/null
+++ b/examples/matlab/old/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}(1:2:end), raw_data.data{p}(2:2:end)), 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
--
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