% ATHOS a Matlab program for the numerical simulation 
% of magnetic suspetibility artifacts in magnetic
% resonance imaging (MRI) for a hollow cylinder geometry.
% 
% (c) Stephane Balac - MAPLY - INSA Lyon - mars 2001
clear; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESCRIPTION OF THE EXPERIMENT
% -  the following variables may be changed
%
global r1;  r1 = .5e-2;            % radius of the 2 cylinders in meter
global r2;  r2 = 1e-2;    
global Xm;  Xm = [1e-4, -1e-5, 0]; % magnetic susceptibilities in the 3 domains
global rho; rho = [0, 10, 0];       % proton density in the 3 domains 
global B0;  B0 = 1;                % strenght of the field B_0 in tesla
global ec;  ec = 3e-3;             % width of the slice in meter
global  g;  g = [1e-2 1e-2 1e-2];  % gradients strength in T/m
beta = 0;                          % angle between the x axis and the
                                   % read-out gradient direction
		          % length of the slice: 
global Lr; Lr = 6e-2;     % along the read-out gradient direction in meter
global Lp; Lp = 6e-2;     % along the phase gradient direction in meter
Np = 64;                % number of pixels = Np^2
%
% -  DO NOT modify the following variables
%
global adim; adim = 1E-2;  % carateristic distance  
global    C;  C = [0 0 0]; % origin of the slice
global    c;  c = sqrt(sum(C.^2)); 
global   ns; ns = [0 0 1];                 % direction of the slice gradient
global   nr; nr = [cos(beta) sin(beta) 0]; % direction of the read-out gradient
global np; np = [-sin(beta) cos(beta) 0];  % direction of the phase gradient
global mu; mu = 1+Xm;
global mu0; mu0 = 4*pi*1.E-7;  % magnetic permeability of vaccum
H0  = B0/mu0;                  % magnetic field strenght
Z = [1 0 0];                   % Direction of B0 in the laboratory frame
% constants for the computation of the magnetic field induced by the cylinder 
global CD; global CC; global CB; global CA; 
if (mu(2)~=mu(1)) % 2 cylinders with distincts magnetic susceptibilities
  CD = 2*r1^2*r2^2*H0*(mu(2)-mu(1))/...
(-r2^2*(mu(1)+mu(2)+mu(1)*mu(2)+mu(2)^2)+r1^2*(mu(2)^2+mu(1)-mu(2)-mu(1)*mu(2)));
  CC = H0 +CD*(mu(2)+mu(1))/((mu(2)-mu(1))*r1^2);
  CA = CC+CD/r1^2;
  CB = CC*r2^2+CD;
else % only one medium
   CA = (mu(1) - 1)*H0/(mu(1)+1);
   CB = (mu(1) - 1)*H0/(r2^2*(mu(1)+1));
   CC = CA;
   CD = 0;
end
% quadrature formula for integration over the slice width
NG = 5;
Gpts = [.906179, .538469, 0., -0.538469, -0.906179];
Gpds = [.236926, .478628, .568888, .478628, .236926];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%warning off % suppresses all subsequent warning messages
times1=cputime;
indi = 0; inde = 0;
IMG=[];       % array containing the pixels intensity
%
% Loop over the pixels
%
for ir = 1:Np
  disp([' ** ',num2str(round((100*ir)/Np)),'% '])
  Pix(1) = - (Lr/2) + (ir-1)*Lr/(Np-1);
for ip = 1:Np
  Pix(2) = - (Lp/2) + (ip-1)*Lp/(Np-1);
  IMG(ir,ip)=0;

for is = 1:NG % Loop over the quadrature nodes
   Pix(3) = ec*Gpts(is)/2;

% Search of the voxels connected to the pixel Pix
   Vox = voxel(Pix(1),Pix(2),Pix(3));
   nvox = size(Vox,2); 
   for ivox=1:nvox
% Computation of the magnetic field induced by the ball in the point Vox(ivox)
     BdB = calBdB(Vox(:,ivox));
     Jac=1./abs(1+BdB(2)/g(1));
     NpA = 40;
     IMG(ir,ip) = IMG(ir,ip) + Gpds(is)*Jac*BdB(3);
     if ((is == 3) & (mod(ir,round(Np/NpA))==0) & (mod(ip,round(Np/NpA))==0))
       inde = inde+1;
       DEF(:,inde) = Vox(:,ivox);
     end 
     if ((is == 3) & (mod(ir,round(Np/NpA))==0) & (mod(ip,round(Np/NpA))==0))
       indi = indi+1;
       DIS(:,indi) = [Pix(1)+nr(1)*BdB(1)/g(1); Pix(2)+nr(2)*BdB(1)/g(1)];
     end    
   end
end
end
end
times2=cputime-times1; disp(['CPU time =',num2str(times2),'sec.']);
%
% Image display
%
hx = Lr/(Np-1);
hy = Lp/(Np-1);
X  = -0.5*Lr:hx:0.5*Lr;
Y  = -0.5*Lp:hy:0.5*Lp;
graylevel = 0.01; % Gray level
Mn = min(min(IMG));
Mx = max(max(IMG));
XX = graylevel*((Mx+Mn)/2)*ones(size(IMG));
IM = min(IMG,XX);
colormap(gray)
imagesc(X,Y,IMG)
title('ATHOS : simulation of the susceptibility artefact')
text(0, 0.6*Lr, '(c 2001) S. Balac - MAPLY - INSA de Lyon')         
axis equal

figure
plot3(DEF(1,:),DEF(2,:),DEF(3,:),'.')
grid
axis equal
title('ATHOS : slice distorsion')
xlabel('read-out gradient direction')
ylabel('phase gradient direction')
zlabel('slice distorstion in meters')

figure
plot(DIS(1,:),DIS(2,:),'.')
grid
axis equal
title('ATHOS : geometric distorsion')
xlabel('read-out gradient direction')
ylabel('phase gradient direction')


