Tutorial: new MRI Loop Coil Tutorial
will need an update to use a human body model Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>pull/10/head
Thorsten Liebig
parent
3a66a02550
commit
3ef6d8df7e
|
|
@ -0,0 +1,296 @@
|
|||
%
|
||||
% Tutorials / 7T MRI Loop Coil
|
||||
%
|
||||
% Describtion at:
|
||||
% http://openems.de/index.php/Tutorial:_MRI_Loop_Coil
|
||||
%
|
||||
% Tested with
|
||||
% - openEMS v0.0.31
|
||||
%
|
||||
% (C) 2013 Thorsten Liebig <thorsten.liebig@uni-due.de>
|
||||
|
||||
close all
|
||||
clear
|
||||
clc
|
||||
|
||||
%% setup the simulation
|
||||
physical_constants; %get some physical constants like c0 and MUE0
|
||||
unit = 1e-3; % all length in mm
|
||||
|
||||
% Loop-Coil parameter
|
||||
loop.length = 80; % length of the loop (in z-direction)
|
||||
loop.width = 60; % width of the loop (in y-direction)
|
||||
loop.strip_width = 5; % metal strip width
|
||||
loop.strip_N_cells = 3; % number of cells over the strip length
|
||||
loop.air_gap = loop.strip_width/3; % air gap width for lumped capacitors
|
||||
loop.pos_x = 110; % position of loop
|
||||
loop.C_gap = 5e-12; % lumped cap value
|
||||
loop.port_R = 10; % feeding port resistance
|
||||
|
||||
%% define the phantom material stackup material
|
||||
materials.name{1}='skin';
|
||||
materials.rad(1)=100;
|
||||
materials.eps_r(1)=49.8;
|
||||
materials.kappa(1)=0.64;
|
||||
materials.density(1)=1100;
|
||||
materials.prio(1)=10;
|
||||
|
||||
materials.name{2}='headbone';
|
||||
materials.rad(2)=95;
|
||||
materials.eps_r(2)=13.4;
|
||||
materials.kappa(2)=0.08;
|
||||
materials.density(2)=1990;
|
||||
materials.prio(2)=11;
|
||||
|
||||
materials.name{3}='CSF';
|
||||
materials.rad(3)=90;
|
||||
materials.eps_r(3)=72.734;
|
||||
materials.kappa(3)=2.2245;
|
||||
materials.density(3)=1007;
|
||||
materials.prio(3)=12;
|
||||
|
||||
materials.name{4}='brain'; % average of white/grey matter
|
||||
materials.rad(4)=86;
|
||||
materials.eps_r(4)=60;
|
||||
materials.kappa(4)=0.69;
|
||||
materials.density(4)=1039;
|
||||
materials.prio(4)=13;
|
||||
|
||||
%% some mesh parameter
|
||||
mat_mesh = 2; % mesh inside the phantom (2mm)
|
||||
Air_Box = 150; % size of the surrounding air box (150mm)
|
||||
|
||||
%% setup FDTD parameter & excitation function
|
||||
% init FDTD structure
|
||||
FDTD = InitFDTD( 'EndCriteria', 1e-4 );
|
||||
|
||||
% define gaussian pulse excitation signal
|
||||
f0 = 300e6; % center frequency
|
||||
fc = 300e6; % 20 dB corner frequency
|
||||
FDTD = SetGaussExcite( FDTD, f0, fc );
|
||||
|
||||
% setup boundary conditions
|
||||
BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
|
||||
FDTD = SetBoundaryCond( FDTD, BC );
|
||||
|
||||
%% setup CSXCAD geometry & mesh
|
||||
CSX = InitCSX();
|
||||
|
||||
%% create loop
|
||||
% setup all properties needed
|
||||
CSX = AddMetal( CSX, 'loop' );
|
||||
CSX = AddLumpedElement( CSX, 'caps_y', 1, 'C', loop.C_gap);
|
||||
CSX = AddLumpedElement( CSX, 'caps_z', 2, 'C', loop.C_gap);
|
||||
|
||||
% horizontal (y-direction) strips
|
||||
start = [loop.pos_x -loop.width/2 -loop.length/2];
|
||||
stop = [loop.pos_x -loop.air_gap/2 -loop.length/2+loop.strip_width];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x -loop.width/2 loop.length/2 ];
|
||||
stop = [loop.pos_x -loop.air_gap/2 loop.length/2-loop.strip_width];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x loop.width/2 -loop.length/2];
|
||||
stop = [loop.pos_x loop.air_gap/2 -loop.length/2+loop.strip_width];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x loop.width/2 loop.length/2 ];
|
||||
stop = [loop.pos_x loop.air_gap/2 loop.length/2-loop.strip_width];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
% vertical (z-direction) strips
|
||||
start = [loop.pos_x -loop.width/2 -loop.length/2+loop.strip_width];
|
||||
stop = [loop.pos_x -loop.width/2+loop.strip_width -loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x -loop.width/2 loop.length/2-loop.strip_width];
|
||||
stop = [loop.pos_x -loop.width/2+loop.strip_width loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x loop.width/2 -loop.length/2+loop.strip_width];
|
||||
stop = [loop.pos_x loop.width/2-loop.strip_width -loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
start = [loop.pos_x loop.width/2 loop.length/2-loop.strip_width ];
|
||||
stop = [loop.pos_x loop.width/2-loop.strip_width loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'loop',10,start,stop);
|
||||
|
||||
% add the lumped capacities
|
||||
start = [loop.pos_x -loop.width/2+loop.strip_width/2-loop.air_gap/2 -loop.air_gap/2];
|
||||
stop = [loop.pos_x -loop.width/2+loop.strip_width/2+loop.air_gap/2 +loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'caps_z',10,start,stop);
|
||||
|
||||
start = [loop.pos_x loop.width/2-loop.strip_width/2-loop.air_gap/2 -loop.air_gap/2];
|
||||
stop = [loop.pos_x loop.width/2-loop.strip_width/2+loop.air_gap/2 +loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'caps_z',10,start,stop);
|
||||
|
||||
start = [loop.pos_x -loop.air_gap/2 loop.length/2-loop.strip_width/2-loop.air_gap/2];
|
||||
stop = [loop.pos_x +loop.air_gap/2 loop.length/2-loop.strip_width/2+loop.air_gap/2];
|
||||
CSX = AddBox(CSX,'caps_y',10,start,stop);
|
||||
|
||||
% add a lumped port as excitation
|
||||
start = [loop.pos_x -loop.air_gap/2 -loop.length/2+loop.strip_width/2-loop.air_gap/2];
|
||||
stop = [loop.pos_x +loop.air_gap/2 -loop.length/2+loop.strip_width/2+loop.air_gap/2];
|
||||
[CSX port] = AddLumpedPort(CSX, 100, 1, loop.port_R, start, stop, [0 1 0], true);
|
||||
|
||||
%% define materials in a loop
|
||||
for n=1:numel(materials.name)
|
||||
CSX = AddMaterial( CSX, materials.name{n} );
|
||||
CSX = SetMaterialProperty( CSX, materials.name{n}, 'Epsilon', materials.eps_r(n), 'Kappa', materials.kappa(n), 'Density', materials.density(n));
|
||||
CSX = AddSphere( CSX, materials.name{n}, materials.prio(n), [0 0 0], materials.rad(n),'Transform',{'Scale',[1 0.8 1] } );
|
||||
end
|
||||
|
||||
%% finalize mesh
|
||||
% create loop mesh (detect only metal and the lumped elements)
|
||||
mesh = DetectEdges(CSX, [], 'SetPropertyType', {'Metal','LumpedElement'});
|
||||
|
||||
% add a dense homegeneous mesh inside the phantom
|
||||
mesh.x = [mesh.x -materials.rad -mat_mesh/2 mat_mesh/2 materials.rad];
|
||||
mesh.y = [mesh.y -materials.rad*0.8 materials.rad*0.8];
|
||||
mesh.z = [mesh.z -materials.rad materials.rad];
|
||||
|
||||
% smooth the mesh for the loop & phantom
|
||||
mesh = SmoothMesh(mesh, mat_mesh);
|
||||
|
||||
% add air spacer
|
||||
mesh.x = [-Air_Box+mesh.x(1) mesh.x mesh.x(end)+Air_Box];
|
||||
mesh.y = [-Air_Box+mesh.y(1) mesh.y mesh.y(end)+Air_Box];
|
||||
mesh.z = [-Air_Box+mesh.z(1) mesh.z mesh.z(end)+Air_Box];
|
||||
|
||||
mesh = SmoothMesh(mesh, c0 / (f0+fc) / unit / 10, 1.5, 'algorithm', 1);
|
||||
|
||||
%% Add Dump boxes (3D box) for E,J and SAR
|
||||
CSX = AddDump(CSX,'Hf','DumpType',11,'FileType',1,'Frequency',f0);
|
||||
CSX = AddDump(CSX,'SAR','DumpType',20,'DumpMode',2,'FileType',1,'Frequency',f0);
|
||||
start = [-120 -120 -120];
|
||||
stop = [+120 +120 +120];
|
||||
CSX = AddBox(CSX,'Hf',0,start,stop);
|
||||
CSX = AddBox(CSX,'SAR',0,start,stop);
|
||||
|
||||
%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
|
||||
start = [mesh.x(1) mesh.y(1) mesh.z(1)];
|
||||
stop = [mesh.x(end) mesh.y(end) mesh.z(end)];
|
||||
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop,'Frequency',f0,'OptResolution',c0/f0/unit/20);
|
||||
|
||||
%% add 10 lines in all direction to make space for PML or MUR absorbing
|
||||
%% boundary conditions
|
||||
mesh = AddPML(mesh, 10);
|
||||
|
||||
%% finaly define the FDTD mesh grid
|
||||
disp(['number of cells: ' num2str(1e-6*numel(mesh.x)*numel(mesh.y)*numel(mesh.z)) ' Mcells'])
|
||||
CSX = DefineRectGrid( CSX, unit, mesh );
|
||||
|
||||
%% prepare simulation folder
|
||||
Sim_Path = mfilename;
|
||||
Sim_CSX = [Sim_Path '.xml'];
|
||||
|
||||
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
|
||||
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
|
||||
|
||||
%% write openEMS compatible xml-file
|
||||
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
|
||||
|
||||
%% show the structure and export as vtk data automatically
|
||||
CSXGeomPlot( [Sim_Path '/' Sim_CSX] , ['--export-polydata-vtk=' Sim_Path]);
|
||||
|
||||
%% run openEMS
|
||||
RunOpenEMS( Sim_Path, Sim_CSX);
|
||||
|
||||
%% postprocessing & do the plots
|
||||
freq = linspace( f0-fc, f0+fc, 501 );
|
||||
port = calcPort(port, Sim_Path, freq);
|
||||
|
||||
Zin = port.uf.tot ./ port.if.tot;
|
||||
s11 = port.uf.ref ./ port.uf.inc;
|
||||
|
||||
P_in = real(0.5 * port.uf.tot .* conj(port.if.tot)); % antenna feed power
|
||||
% get the feeding power for frequency f0
|
||||
P0_in = interp1(freq, P_in, f0);
|
||||
|
||||
%%
|
||||
% plot reflection coefficient S11
|
||||
figure
|
||||
h = plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
|
||||
grid on
|
||||
title( 'reflection coefficient S_{11}' );
|
||||
xlabel( 'frequency f / MHz' );
|
||||
ylabel( 'reflection coefficient |S_{11}| (dB)' );
|
||||
|
||||
% plot feed point admittance
|
||||
figure
|
||||
h = plot( freq/1e6, real(1./Zin), 'k-', 'Linewidth', 2 );
|
||||
hold on
|
||||
grid on
|
||||
plot( freq/1e6, imag(1./Zin), 'r--', 'Linewidth', 2 );
|
||||
title( 'feed port admittance' );
|
||||
xlabel( 'frequency f (MHz)' );
|
||||
ylabel( 'admittance Y_{in} (S)' );
|
||||
legend( 'real', 'imag' );
|
||||
|
||||
%% read SAR values on a xy-plane (range)
|
||||
[SAR SAR_mesh] = ReadHDF5Dump([Sim_Path '/SAR.h5'], 'Range', {[],[],0});
|
||||
SAR = SAR.FD.values{1};
|
||||
|
||||
%%
|
||||
% SAR plot
|
||||
figure()
|
||||
[X Y] = ndgrid(SAR_mesh.lines{1},SAR_mesh.lines{2});
|
||||
colormap('hot');
|
||||
h = pcolor(X,Y,(squeeze(SAR)));
|
||||
% h = pcolor(X,Y,log10(squeeze(SAR)));
|
||||
set(h,'EdgeColor','none');
|
||||
xlabel('x -->');
|
||||
ylabel('y -->');
|
||||
title('local SAR');
|
||||
axis equal tight
|
||||
|
||||
%%
|
||||
[H_field H_mesh] = ReadHDF5Dump([Sim_Path '/Hf.h5'], 'Range',{[],[],0});
|
||||
% calc Bx,By, B1p, B1m normalize to the input-power
|
||||
Bx = MUE0*H_field.FD.values{1}(:,:,1,1)/sqrt(P0_in);
|
||||
By = MUE0*H_field.FD.values{1}(:,:,1,2)/sqrt(P0_in);
|
||||
B1p = 0.5*(Bx+1j*By);
|
||||
B1m = 0.5*(Bx-1j*By);
|
||||
% create a 2D grid to plot on
|
||||
[X Y] = ndgrid(H_mesh.lines{1},H_mesh.lines{2});
|
||||
|
||||
filter = sqrt(X.^2+Y.^2)<0.1;
|
||||
Dump2VTK([Sim_Path '/B1p_xy.vtk'], abs(B1p).*filter, H_mesh, 'B-Field');
|
||||
Dump2VTK([Sim_Path '/B1m_xy.vtk'], abs(B1m).*filter, H_mesh, 'B-Field');
|
||||
|
||||
% B1+ plot
|
||||
figure()
|
||||
h = pcolor(X,Y,log10(abs(B1p)));
|
||||
set(h,'EdgeColor','none');
|
||||
xlabel('x -->');
|
||||
ylabel('y -->');
|
||||
title('B_1^+ field (dB)');
|
||||
axis equal tight
|
||||
|
||||
% B1- plot
|
||||
figure()
|
||||
h = pcolor(X,Y,log10(abs(B1m)));
|
||||
set(h,'EdgeColor','none');
|
||||
xlabel('x -->');
|
||||
ylabel('y -->');
|
||||
title('B_1^- field (dB)');
|
||||
axis equal tight
|
||||
|
||||
%% dump to vtk to view in Paraview
|
||||
Dump2VTK([Sim_Path '/SAR_xy.vtk'], SAR, SAR_mesh, 'SAR');
|
||||
|
||||
%%
|
||||
ConvertHDF5_VTK([Sim_Path '/Hf.h5'],[Sim_Path '/B1_xy'], 'Range',{[],[],0}, 'weight', MUE0/sqrt(P0_in), 'FieldName', 'B1-field');
|
||||
|
||||
%% validation by calculating the power budget
|
||||
% calculate the radiated power
|
||||
nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, 0, 0);
|
||||
p_rad = nf2ff.Prad;
|
||||
|
||||
% read dissipated power in the phantom
|
||||
p_loss = ReadHDF5Attribute([Sim_Path '/SAR.h5'],'/FieldData/FD/f0','power');
|
||||
|
||||
% display results, should add up to 100% (+/- 5% error margin)
|
||||
disp([' power loss in phantom: ' num2str(p_loss) ' W (' num2str(p_loss/P0_in*100) '%)']);
|
||||
disp([' power radiated: ' num2str(p_rad) ' W (' num2str(p_rad/P0_in*100) '%)']);
|
||||
Loading…
Reference in New Issue