-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathdemo.m
186 lines (156 loc) · 5.13 KB
/
demo.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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
clear all; close all; clc;
addpath('utils/')
addpath(genpath('fwtoolbox_v1_code/')) % ISMRM water-fat toolbox
%% load data
% load mat-file with imaging data after reconstruction
load('20151101_151725_0302_ImDataParams.mat')
signal0 = ImDataParams.signal;
signal = signal0;
TE_s = ImDataParams.TE_s;
centerFreq_Hz = ImDataParams.centerFreq_Hz;
B0dir = ImDataParams.B0dir;
voxelSize_mm = ImDataParams.voxelSize_mm;
isl = 37;
p0 = angle(signal(:, :, isl, end));
% load mat-file with coefficients of the spherical harmonic expantion of
% the magnet inhomogeneities and the shimfield
% Sorry, I am not allowed to add this file to the repository
load('20151101_151725_0302_B0params.mat')
B0params
%% estimate and demodulate field contributions
% magnet inhomogeneities
sz = size(signal);
matrixSize = sz(1:3);
transform = B0params.affMat_ijk2xyMz;
magnetInhomogeneities_Hz = get_magnetInhomogeneities_Hz(B0params, matrixSize, transform);
signal = demodulate_field_Hz(signal, TE_s, magnetInhomogeneities_Hz);
p1 = angle(signal(:, :, isl, end));
% shim field
shimField_Hz = get_shimField_Hz(B0params.shimValues, matrixSize, transform);
signal = demodulate_field_Hz(signal, TE_s, shimField_Hz);
p2 = angle(signal(:, :, isl, end));
% object-based fieldmap estimate
objectBasedFieldmap_Hz = get_objectBasedFieldmap_Hz(signal, voxelSize_mm, B0dir, centerFreq_Hz);
signal = demodulate_field_Hz(signal, TE_s, objectBasedFieldmap_Hz);
p3 = angle(signal(:, :, isl, end));
% residual linear field
cropDim = 1;
residualLinearField_Hz = get_residualLinearField_Hz(signal, TE_s, cropDim);
signal = demodulate_field_Hz(signal, TE_s, residualLinearField_Hz);
p4 = angle(signal(:, :, isl, end));
%% water--fat separation
imDataParams.images = reshape(signal(:, :, isl, :), ...
[sz(1), sz(2), 1, 1, sz(4)]); % add dummy coil dimension, TE dimension last
imDataParams.TE = TE_s;
imDataParams.FieldStrength = 3;
imDataParams.PrecessionIsClockwise = 1; % Philips scanner follows convention of clockwise precession
algoParams.species(1).name = 'water';
algoParams.species(1).frequency = 0;
algoParams.species(1).relAmps = 1;
algoParams.species(2).name = 'fat (7 peaks)';
algoParams.species(2).frequency = [-3.8, -3.4, -3.1, -2.68, -2.46, -1.95, -0.5, 0.49, 0.59];
algoParams.species(2).relAmps = [0.0899, 0.5834, 0.0599, 0.0849, 0.0599, 0.0150, 0.0400, 0.01, 0.0569];
% hIDEAL
algoParams.Visualize = 0;
algoParams.AlwaysShowGUI = 0;
outParams_hIDEAL = fw_i2cm0c_3pluspoint_tsaojiang(imDataParams, algoParams);
% graph cut
algoParams.range_fm = [-1000, 1000];% Range of field map values
algoParams.NUM_FMS = 500;% Number of field map values to discretize
outParams_graphcut = fw_i2cm1i_3pluspoint_hernando_graphcut(imDataParams, algoParams);
% without demodulation steps for comparison
imDataParams.images = reshape(signal0(:, :, isl, :), ...
[sz(1), sz(2), 1, 1, sz(4)]); % add dummy coil dimension, TE dimension last
outParams0_hIDEAL = fw_i2cm0c_3pluspoint_tsaojiang(imDataParams, algoParams);
outParams0_graphcut = fw_i2cm1i_3pluspoint_hernando_graphcut(imDataParams, algoParams);
%% plot results
% standard WFI results withouth any demodulation steps
figure('position', [0, 0, 1000, 1000])
colormap gray
subplot(2, 2, 1)
imagesc(abs(outParams0_hIDEAL.species(1).amps))
set(gca, 'xtick', [])
set(gca, 'ytick', [])
ylabel('water')
title('hIDEAL')
subplot(2, 2, 2)
imagesc(abs(outParams0_graphcut.species(1).amps))
axis off
title('graph cut')
subplot(2, 2, 3)
imagesc(abs(outParams0_hIDEAL.species(2).amps))
set(gca, 'xtick', [])
set(gca, 'ytick', [])
ylabel('fat')
subplot(2, 2, 4)
imagesc(abs(outParams0_graphcut.species(2).amps))
axis off
saveas(gcf, 'stdWFI.png')
% WFI results after proposed demodulation steps
figure('position', [1000, 0, 1000, 1000])
colormap gray
subplot(2, 2, 1)
imagesc(abs(outParams_hIDEAL.species(1).amps))
set(gca, 'xtick', [])
set(gca, 'ytick', [])
ylabel('water')
title('hIDEAL')
subplot(2, 2, 2)
imagesc(abs(outParams_graphcut.species(1).amps))
axis off
title('graph cut')
subplot(2, 2, 3)
imagesc(abs(outParams_hIDEAL.species(2).amps))
set(gca, 'xtick', [])
set(gca, 'ytick', [])
ylabel('fat')
subplot(2, 2, 4)
imagesc(abs(outParams_graphcut.species(2).amps))
axis off
saveas(gcf, 'proposedWFI.png')
% row 1: magnitude image TE1, estimated field contributions in [Hz]
% row 2: phase images of TE3 after demodulation of field contribution in row 1
figure('position', [0, 0, 2500, 1000])
colormap('gray')
subplot(2, 5, 1)
imagesc(abs(signal(:, :, isl, 1)))
axis off
title('magnitude TE1')
subplot(2, 5, 2)
imagesc(magnetInhomogeneities_Hz(:, :, isl))
colorbar
axis off
title('B0 inhom.')
subplot(2, 5, 3)
imagesc(shimField_Hz(:, :, isl))
title('shim')
colorbar
axis off
subplot(2, 5, 4)
imagesc(objectBasedFieldmap_Hz(:, :, isl))
title('OBFFME')
colorbar
axis off
subplot(2, 5, 5)
imagesc(residualLinearField_Hz(:, :, isl))
title('res. lin. field.')
colorbar
axis off
subplot(2, 5, 6)
imagesc(p0)
ylabel('phase TE3')
set(gca, 'xtick', [])
set(gca, 'ytick', [])
subplot(2, 5, 7)
imagesc(p1)
axis off
subplot(2, 5, 8)
imagesc(p2)
axis off
subplot(2, 5, 9)
imagesc(p3)
axis off
subplot(2, 5, 10)
imagesc(p4)
axis off
saveas(gcf, 'field_contributions.png')