-
Notifications
You must be signed in to change notification settings - Fork 2
/
simualteSignalExample.m
112 lines (83 loc) · 3.62 KB
/
simualteSignalExample.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
% This example simulates field perturbation and corresponding multi GRE
% signal from WM models
% Two single models are provided, one 2D model (by defaut) and one 3D model
% (that you can uncomment to run)
% You can load your own WM model create by createOneWMModelExample.m
clear
% close all
your_folder = '/home/rhedouin/Code/Whist';
% location of the toolbox
addpath(genpath(your_folder))
%%%%%%%%%%%% Load a WM model with a single 2D axon
model_path = '/home/rhedouin/Code/Whist/data/oneAxon2D.mat';
%%%%%%%%%%%% Load your WM model
% model_path = '/home/rhedouin/Code/Whist/WMmodel/MyWMModel.mat';
%%%%%%%%%%%% Load a WM model with a single 3D axon
% model_path = '/project/3015069.04/code/Whist/data/oneAxon3D.mat';
load(model_path)
number_dims = ndims(mask);
model_parameters.dims = size(mask);
plot_model = 1;
% Create and plot your model
[model, zoomed_model, FVF, g_ratio] = createModelFromData(axon_collection, mask, plot_model);
%%%%%%%%%%% Set parameters
% mask (required)
model_parameters.mask = mask;
% myelin (required: T2, xi, xa)
model_parameters.myelin.weight= 0.5;
model_parameters.myelin.T2 = 15*1e-3;
model_parameters.myelin.T1 = 500*1e-3;
model_parameters.myelin.proton_density= 0.5;
model_parameters.myelin.xi = -0.1; % myelin anisotropic susceptibility (ppm)
model_parameters.myelin.xa = -0.1; % myelin isotropic susceptibility (ppm)
% intra axonal (required: T2)
model_parameters.intra_axonal.weight = 1;
model_parameters.intra_axonal.T2 = 50*1e-3;
model_parameters.intra_axonal.T1 = 1.5;
model_parameters.intra_axonal.proton_density= 1;
model_parameters.intra_axonal.xi= 0;
% extra axonal (required: T2)
model_parameters.extra_axonal.weight = 1;
model_parameters.extra_axonal.T2 = 50*1e-3;
model_parameters.extra_axonal.T1 = 1.5;
model_parameters.extra_axonal.proton_density= 1;
model_parameters.extra_axonal.xi= 0;
% main magnetic field strength in Tesla (required)
model_parameters.B0 = 3;
% magnetic field orientation (required)
model_parameters.theta = pi/2;
model_parameters.phi = 0;
model_parameters.field_direction = [sin(model_parameters.theta)*cos(model_parameters.phi), ...
sin(model_parameters.theta)*sin(model_parameters.phi), ...
cos(model_parameters.theta)];
% TE (required)
model_parameters.TE = (2:3:80)*1e-3;
% optional, needed to include T1 effect in signal weights
model_parameters.flip_angle = 20;
model_parameters.TR = 60*1e-3;
%
model_parameters.include_proton_density = 1;
model_parameters.include_T1_effect = 1;
% apply or not the Lorentzian correction (see He and Yablonskiy (2009))
model_parameters.Lorentzian = 0;
%%%%%%%%%%%
% Estimate the relative weights of each compartment
% model_parameters = computeCompartmentSignalWeight(model_parameters);
%%%%%%%%%% Simulate the field perturbation from the WM model and the multi GRE signals
if model_parameters.Lorentzian
[signal_original, field] = simulateSignalFromModelWithLorentzianCorrection(axon_collection, model_parameters);
elseif model_parameters.Lorentzian == 0
[signal_original, field] = simulateSignalFromModel(axon_collection, model_parameters);
else
error('Lorentzian parameter need to be 0 or 1');
end
%%%%%%%%%% Plot frequency histogramm, field and signals
options.mask = mask;
createHistogramFieldPerturbation(model, field, options);
if number_dims == 2
plot2DFieldAndSignal(field, signal_original.total, model_parameters.TE, model_parameters.field_direction)
elseif number_dims == 3
plot3DFieldAndSignal(field, signal_original.total, model_parameters.TE, model_parameters.field_direction)
else
error('dimension of the model should be 2 or 3');
end