-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsimpleRFanimations.m
181 lines (147 loc) · 5.25 KB
/
simpleRFanimations.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
%%% simpleRFanimations
%%% ---
%%% Code to create animated visualisations of spin dynamics for a Hahn
%%% Echo, a simple spin echo and a stimulated echo.
%%%
%%% Code was first created years ago. Put on GitHub for easier sharing June
%%% 2019.
%%%
%%% Daniel Gallichan
saveMovie = 0; % set to one to save out animation frames
savedir='~/temp/mov'; % temporary folder to save animation frames to
addpath(genpath('tools'))
%%% Hahn Echo:
hahnopts.Angles = [90 90];
hahnopts.Flip_times = [0, 1];
hahnopts.rot_funcs(1:2) = {'rotx'};
hahnopts.t = linspace(-.05,2,501);
hahnopts.z_max = .5; hahnopts.z_min =-hahnopts.z_max;
hahnopts.flip_images = 40; % no of frames during RF pulse itself
% % Simple Spin Echo:
SEopts.Angles = [90 180];
SEopts.Flip_times = [0, 1];
SEopts.rot_funcs(1) = {'rotx'};
SEopts.rot_funcs(2) = {'roty'};
SEopts.t = linspace(-.05,3.5,501);
SEopts.z_max = .5; SEopts.z_min = -SEopts.z_max;
SEopts.noSpins = 500;
SEopts.flip_images = 40;
%%% Stimulated echo
stimEcho.Angles = [90 90 90];
stimEcho.Flip_times = [0, .5, 1.5];
stimEcho.rot_funcs(1:3) = {'rotx'};
stimEcho.t = linspace(-0.05,3.5,501);
stimEcho.z_max = 5; stimEcho.z_min = -stimEcho.z_max;
stimEcho.noSpins = 500;
stimEcho.flip_images = 40;
sim = simEvolution(hahnopts);
% sim = simEvolution(SEopts); % <-- use this line for SE
% sim = simEvolution(stimEcho); % <-- use this line for stimulated echo
M = sim.M; newM = sim.newM; newtime_i = sim.newtime_i; t = sim.t; rf_tracking = sim.rf_tracking; rot_funcs = sim.rot_funcs; Flip_times = sim.Flip_times; Angles = sim.Angles;
meanSignalY = mean(M(2,:,:),3);
meanSignalX = mean(M(1,:,:),3);
x0 = [-1 0 0; 1 0 0];
y0 = [ 0 -1 0; 0 1 0];
z0 = [ 0 0 -1; 0 0 1];
figure
set(gcf,'color','w')
set(gcf,'Position',[ 241 98 758 840])
TopAxis = subplot(211);
BottomAxis = subplot(212);
set(TopAxis,'Position',[ -0.068 0.005 1.17 1.17]);
set(BottomAxis,'Position',[ 0.0354 0.0435 0.9511 0.1714 ]);
subplot(TopAxis)
Az = +37.5-90; El = 30;
view(Az,El);
lw = 1.5;
M = newM;
imCount = 1;
rf_flash = 1;
no_flips = length(Flip_times);
animate = 1;
animSkip = 4;
if animate
for dt = 1:animSkip:size(M,2)
subplot(TopAxis)
cla
x = squeeze(M(1,dt,:));
y = squeeze(M(2,dt,:));
z = squeeze(M(3,dt,:));
vertices = [0 0 0; x y z];
faces = ones(length(x)-1,3);
faces(:,2) = 2:length(x);
faces(:,3) = 3:length(x)+1;
colordata = repmat(linspace(0,1,length(faces))',[1 3]);
colordata(:,1) = colordata(end:-1:1,1);
colordata(:,2) = 0;
patch('vertices',vertices,'faces',faces,'facevertexcdata',colordata,'facecolor','flat','edgecolor','none')
view(Az,El)
axSize = 1.2;
showAxis3d(axSize)
set(gca,'DataAspectRatioMode','Manual','DataAspectRatio',[1 1 1])
axis vis3d
axis off
hold all
linesToPlot = vertices(2:10:end,:);
nLP = size(linesToPlot,1);
zLP = zeros(nLP,1);
for iL = 1:nLP
line([0 linesToPlot(iL,1)],[0 linesToPlot(iL,2)],[0 linesToPlot(iL,3)],'linewidth',1,'color','k')
end
rf_flash = -rf_flash;
if rf_flash > 0
aColor = [13 151 21]/255;
rotlstyle = {'linewidth',5,'color',aColor};
else
aColor = [13 90 21]/255;
rotlstyle = {'linewidth',5,'color',aColor};
end
if (rf_flash) && (rf_tracking(dt))
switch rf_tracking(dt)
case 1
aFinish = [1 0 0];
aNormal = [0 0 1];
case 2
aFinish = [0 1 0];
aNormal = [0 0 1];
case 3
aFinish = [0 0 1];
aNormal = [1 0 0];
end
arrow([0 0 0],aFinish,50,'BaseAngle',90,'tipangle',30,'width',20,'normaldir',aNormal,'facecolor',aColor)
end
subplot(BottomAxis)
cla
hold on
box on
grid on
set(gca,'linewidth',1)
plot(t,meanSignalY,'linewidth',lw)
plot(t(newtime_i(dt)),meanSignalY(newtime_i(dt)),'o','linewidth',lw)
for flip_i = 1:no_flips
plot([Flip_times(flip_i) Flip_times(flip_i)],[-1 1],'--k','linewidth',2)
rotaxis = char(rot_funcs(flip_i)); rotaxis = rotaxis(end);
text(Flip_times(flip_i), 1.25, [num2str(Angles(flip_i)),'^\circ_', rotaxis],'HorizontalAlignment','center','fontsize',16)
end
if (rf_tracking(dt) > 0) && (rf_flash)
line([t(newtime_i(dt)) t(newtime_i(dt))],[-1 1],rotlstyle{:});
end
line([min(t) max(t)],[0 0],'color','k','linewidth',2)
axis([min(t) max(t) -1.01 1.01])
set(gca,'xticklabel',[],'yticklabel',[])
drawnow
if saveMovie
numString = num2str(imCount);
if length(numString)==1
numString = ['00',numString];
end
if length(numString)==2
numString = ['0',numString];
end
FileName = ['mov',numString];
export_fig([savedir '/' FileName],'-nocrop')
imCount = imCount + 1;
end
end
end