forked from MIITT-MRI-Jianglab/Abdominal_MR_Phantom
-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_phantom_generation.m
164 lines (137 loc) · 5.81 KB
/
main_phantom_generation.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
%% main_phantom_generation
%
% Preparation for Abdominal_MR_Phantom generation:
% (1) user-defined file name
% (2) sequence type
% (3) sampling trajectory
% (4) FOV and spatial resolution
% (5) sampling lines
% (6) respratory motion file and temporal resolution
% (7) coil sensitivity file and coil numbers
% (8) Fat-water chemical shift frequency
%
%%%%%%%%%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define tissue properties and respiratory motion.
% Generate voxelized volumetric mask with indexes.
% ------------------------------------------------
%
%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define sequence parameters.
% Generate signal evolution using Bloch simulation.
% Generate k-space data for each time frame.
% ------------------------------------------------
%
%%%%%%%%%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data reconstruction and gridding
% ------------------------------------------------
%
% -----------------------------------------------------------------------------------------
% Realistic 4D abdominal phantom for magnetic resonance imaging
% Wei-Ching (Tina) Lo
% wxl317@case.edu
% Case Western Reserve University
% July 2018
% -----------------------------------------------------------------------------------------
clear; close all;
% Define file name
savename = 'voxim15';
% Define sequence parameters:
% 'SpoiledGradientEcho' => T1-weighted image
% 'SpoiledGradientEchoWithFatSat' => Perfusion-weighted image
% 'InversionRecoveryLookLocker' => T1 mapping
% 'SingleSpinEcho' => T2 mapping
% 'SingleSpinEchoWithFatSat' => Diffusion-weighted image
% 'MultiEchoSpoiledGradientEcho' => Proton density fat fraction
sigtype = 'SingleSpinEcho';
% Define sampling trajectory:
% 'cartesian'
% 'radial'
% 'spiral'
samptraj = 'cartesian';
% Define 3D resolution and FOV
fov = 420; % field-of-view (mm)
mtx = 256; % matrix size
npar = 64; % # of partitions
slthick = 3; % slice thickness (mm)
nset = 1; % # of sets
% Define number of sampling lines
np_cartesian = 256; % # of Cartesian lines
np_radial = 288; % # of projections/spokes
np_spiral = 48; % # of spiral arms
sampmode = 'demo'; % sampling mode: 'demo', 'simple', 'eachEcho'
% Note that "sampmode" will affect the simulation time (mins to hours, depending on the number of Echoes).
% Use 'eachEcho' option only when Echo-by-Echo simulation is nessasary
% Define 3D respiratory motion curve and temporal resolution (eq. each TR or each volume)
respmotion = 'respmov.mat';
tempres = 400; % temporal resolution (ms)
tempdur = 4000; % duration of each respiratory cycle (ms)
SImov = 13; % largest superior-inferior (SI) excursion (mm)
APmov = 6.5; % largest anterior-posterior (AP) excursion (mm)
LRmov = 2; % largest left-right (LR) excursion (mm)
% Define coil sensitivity map
coilmap = 'origcmap.mat';
nc = 20; % # of coils
% Define fat-water chemical shift
FWshift = 220; % water and fat are separated by approximately 440Hz in a 3T static field
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
datapath = fileparts(mfilename('fullpath'));
cd(datapath)
addpath(genpath(datapath));
% Load defined tissue properties: T1, T2, ADC, PDFF
tissueprop = tissueproperty;
% respiratory motion pattern
load(respmotion); % respmov
tframe = floor(tempdur/tempres);
linv = genresp(respmov,tframe,[SImov APmov LRmov]);
% Convert mesh models to voxels (This step will take mins to hours, depending on the number of time frames)
xpts = -round((fov/420)*mtx/2)+1:round((fov/420)*mtx/2);
ypts = -round((fov/420)*mtx/2)+1:round((fov/420)*mtx/2); % x matrix size = y matrix size
zpts = linspace(round(-npar*slthick/3+35),round(npar*slthick/3+35),npar);
phanimg = mesh2model(tissueprop,linv,xpts,ypts,zpts);
% Show phantom images
showimg(phanimg);title('Indexed phantom images: axial plane');
save([savename '_phanimg.mat'],'phanimg','-v7.3')
fprintf('Phantom generation done\n');
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare FFT/NUFFT
if strcmp(samptraj,'cartesian')
np = np_cartesian;
elseif strcmp(samptraj,'radial')
np = np_radial;
elseif strcmp(samptraj,'spiral')
np = np_spiral;
end
opts = prepareNUFFT(mtx,np,samptraj,'linear_sorted','fov',fov,'FWshift',FWshift);
% Load in-vivo/simulated coil sensitivity maps
load(coilmap);
cmap = gencmap([mtx mtx npar],nc,origcmap);
% Sequence parameters
[seqparam,defseq] = setseqparam(sigtype,[np npar nset],sampmode);
% Bloch simulation
sigevo = gensigevo(tissueprop,seqparam);
% Convert voxels to phantom images (This step will take mins to hours, depending on the number of time frames)
nt = length(defseq.demosig);
[nr,~] = size(opts.kx);
mixsamp = zeros(nr,np,nc,npar,nt,'single');
for itp = 1:nt
imPall = model2voximg(phanimg(:,:,:,mod(defseq.demosig(itp)-1,tframe)+1),sigevo(defseq.demosig(itp),:,:)); % Ground truth images
if itp == 1
nval = calcnoiselvl(imPall,cmap);
end
mixsamp(:,:,:,:,itp) = voximg2ksp(imPall,cmap,nval,opts); % k-space + noise
end
save([savename '_mixsamp.mat'],'mixsamp','-v7.3')
fprintf('Data acquisition done\n');
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert 4D phantom k-space to images
reconimg = ksp2img(mixsamp,opts,cmap);
fprintf('Data reconstruction done\n');
% Show phantom images
showimg(reconimg(:,:,round(npar/2),:));colormap(gray);title('Reconstructed phantom images: axial plane')
showimg(imrotate(squeeze(reconimg(round(mtx/2),:,:,:)),90));colormap(gray);title('Reconstructed phantom images: frontal plane')