/
simplefft.m
executable file
·92 lines (74 loc) · 3.73 KB
/
simplefft.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
classdef simplefft < handle
% Created by Alexander Fyrdahl <alexander.fyrdahl@gmail.com>
methods (Static)
function image = process(conn,config,metadata,log)
log.info('Reconstructing %s\n',config);
try
log.info('Metadata is %s\n', ismrmrd.xml.serialize(metadata));
encoding_x = metadata.encoding.encodedSpace.matrixSize.x;
encoding_y = metadata.encoding.encodedSpace.matrixSize.y;
recon_x = metadata.encoding.reconSpace.matrixSize.x;
recon_y = metadata.encoding.reconSpace.matrixSize.y;
num_coils = metadata.acquisitionSystemInformation.receiverChannels;
catch
disp("Error deserializing XML header");
end
group = ismrmrd.Acquisition;
acqGroup = cell(1,0);
while true
meas = next(conn);
if isa(meas, 'ismrmrd.Acquisition')
% Ignore non-imaging data
if ~(meas.head.flagIsSet(meas.head.FLAGS.ACQ_IS_NOISE_MEASUREMENT) || ...
meas.head.flagIsSet(meas.head.FLAGS.ACQ_IS_PHASECORR_DATA))
acqGroup{end+1} = meas;
% append(group, meas.head, meas.traj(:), meas.data(:));
end
if meas.head.flagIsSet(meas.head.FLAGS.ACQ_LAST_IN_SLICE)
break
end
else
break
end
end
ksp = cell2mat(permute(cellfun(@(x) x.data, acqGroup, 'UniformOutput', false), [1 3 2]));
% ksp = reshape([group.data{:}],encoding_x,num_coils,encoding_y);
ksp = permute(ksp,[1 3 2]);
ksp = fftshift(ifft(fftshift(ksp,1),[],1),1);
im = ksp;
%{
removing for this simple test since pulseq doesn't update this
info and it is cropping incorrectly.
ind1 = floor((encoding_x - recon_x)/2)+1;
ind2 = floor((encoding_x - recon_x)/2)+recon_x;
im = ksp(ind1:ind2,:,:);
%}
im = fftshift(ifft(fftshift(im,2),[],2),2);
im = sqrt(sum(abs(im).^2,3));
im = im.*(32768./max(im(:)));
im = round(im);
im = int16(im);
image = ismrmrd.Image(im);
% find the center Idx fto set the output meta attributes.
kspace_encode_step_1 = cellfun(@(x) x.head.idx.kspace_encode_step_1, acqGroup);
centerLin = cellfun(@(x) x.head.idx.user(6), acqGroup);
centerIdx = find(kspace_encode_step_1 == centerLin, 1);
% field_of_view is mandatory
image.head.field_of_view = single([metadata.encoding(1).reconSpace.fieldOfView_mm.x ...
metadata.encoding(1).reconSpace.fieldOfView_mm.y ...
metadata.encoding(1).reconSpace.fieldOfView_mm.z]);
% Set ISMRMRD Meta Attributes
meta = struct;
meta.DataRole = 'Image';
meta.ImageProcessingHistory = 'MATLAB';
meta.WindowCenter = uint16(16384);
meta.WindowWidth = uint16(32768);
meta.ImageRowDir = acqGroup{centerIdx}.head.read_dir;
meta.ImageColumnDir = acqGroup{centerIdx}.head.phase_dir;
% set_attribute_string also updates attribute_string_len
image = image.set_attribute_string(ismrmrd.Meta.serialize(meta));
conn.send_image(image);
conn.send_close();
end
end
end