% ARAMIS a Matlab program for the numerical SIMulation 
% of metallic implant induced ARtifActs in magnetic resonance
% imaging (MRI) 
%
% (c) S. Balac  - Univ. de Bretagne Occidentale   - 1999 (version 1.0)
% (c) S. Balac  - ENSSAT - Université de Rennes 1 - 2007 (version 2.0)
% and L. Chupin - INSA de Lyon (laurent.chupin@insa-lyon.fr)
%
%
%    Stephane Balac  
%    ENSSAT - Université de Rennes 1
%    6 Rue de Kerampont - BP 80518 
%    22305 LANNION
%    FRANCE 
%    <http://perso.univ-rennes1.fr/stephane.balac> 
%
%   Laurent Chupin
%   Institut Camille Jordan
%   Bât. L. de Vinci, INSA de Lyon
%   69621 VILLEURBANNE cedex, FRANCE
%   <http://maths.insa-lyon.fr/~chupin> 
%
% The code source can be download at:
%    http://perso.univ-rennes1.fr/stephane.balac/irm/aramis
% or by sending a request to:
%    stephane.balac@univ-rennes1.fr
clear
%=====================================================================
% DESCRIPTION OF THE EXPERIMENT - these variables may be changed
%
% Unless otherwise specified, coordinates refer to the laboratory frame
%
% Description of the ball
% -----------------------
% 
% Cu : Xm = -2E-5  sigma = 5.9E7
% Ti : Xm = 2E-4   sigma = 1.9E6
global r0; r0 = 1e-2;          % its radius in meter
global Xm; Xm = -2E-5;         % its magnetic susceptibility
% global Xm; Xm = 0 % uncomment this line to hidde the susceptibility artifact
global sigma; sigma =   5.9E7;  % its electric conductivity
% global sigma; sigma = 0; % uncomment this line to hidde the eddy current artifact
%
% Description of the image sequence (in the laboratory frame)
% ---------------------------------
T1 = 1e-2; % Relaxation times of the surrounding area  
T2 = 1e-3;
global  C;  C = [0 0 0];  % position of the center of the slice
global B0; B0 = 2.5;      % strenght of the field B_0 in tesla
global B1; B1 = 1E-3;     % strenght of the field B_1 in tesla
TE = 1;                   % Echo time
global eta0; eta0 = 1;    % Strengh of the equilibrium magnetisation
global ec; ec = 3e-3;     % width of the slice in meter
global  g;  g = 1e-2;     % intensity of the gradients in tesla per meter
global ns; ns = [0 1 0];  % direction of the slice selection gradient
global nr; nr = [1 0 0];  % direction of the read-out gradient
global np; np = [0 0 -1];  % direction of the phase gradient
global Lr; Lr = 6e-2;     % length of the slice along nr in meter
global Lp; Lp = 6e-2;     % length of the slice along np in meter
           Np = 128;       % number of pixels (Np^2)
%=====================================================================
% DECLARATION OF INTERNAL CONSTANTS - do not change the assigned values
% warning off % suppresses all subsequent warning messages
mu0 = 4*pi*1E-7;                      % permeability of vaccum
global gamma; gamma = 42.58E6;        % gyromagnetic ratio for proton
global omega; omega = gamma*B0;       % Larmor pulsation
global beta; beta = (Xm*B0*r0^3)/3; 
global  c;  c = sqrt(sum(C.^2)); 
NG=5;  % Quadrature formula for integration over the slice width
Gpts = [.906179, .538469, 0., -0.538469, -0.906179];
Gpds = [.236926, .478628, .568888, .478628, .236926];

k = [0 0 1]; % Direction of B_0 in the laboratory frame 
global Z; Z = [dot(k,nr),dot(k,np),dot(k,ns)]; % and in (C,nr,np,ns)

times1=cputime;
IMG=[];       % array containing the pixels intensity
%
% Loop over the pixels
%
for ir = 1:Np % pixels in the read out gradient direction 
  disp([' ** ',num2str(round((100*ir)/Np)),'% '])
  Pix(1) = - (Lr/2) + (ir-1)*Lr/(Np-1);  
for ip = 1:Np % pixels in the phase gradient direction 
  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;

% Computation of the voxels connected to the pixel Pix
   if Xm ~= 0  % susceptibility effects are taken into account
     Vox = voxel(Pix);   % location of the voxels  connected to the current pixel
     nvox = size(Vox,2); % number of voxels connected to the current pixel
     else % susceptibility effects are not taken into account
     Vox(:,1) = Pix';
     nvox = 1; % one voxel is connected to the current pixel
   end

   for ivox=1:nvox
     P = Vox(:,ivox); % coordinates of the voxel in the frame (C,nr,np,ns)
% Computation of the static magnetic field induced by the ball in P
     BdB = calBdB(P);
     Jac=1./abs(1+BdB(2)/g); % Jacobian
% Computation of the Eddy current induced field in P
     if sigma ~= 0  % RF perturbation are taken into account
       [U,V] = calBRF(P);
       UV = norm(U)+norm(V);
       if UV > B1/40 % the RF perturbation is not to small
         M1 = solvebloch(pi/2,[0 0 eta0]',3,T1,T2,U,V);
         M2 = relax(P,M1,TE/2);
         M3 = solvebloch(pi,M2',3,T1,T2,U,V);
         SRF = sqrt(M3(1)^2+M3(2)^2);
       else % RF perturbation are negligible
         SRF = 1;
       end
     else % RF perturbation are not taken into account
       SRF = 1;
     end
% Image intensity distorsions
     IMG(ir,ip) = IMG(ir,ip) + Gpds(is)*Jac*BdB(3)*SRF;
   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;
c = mean(mean(IMG));
ratio = 10;
Mn = min(min(IMG));
Mx = max(max(IMG));
d1 = (Mx-c);
d2 = (c-Mn);
IM = min(IMG,max(0,c+ratio*d1));
IM = max(IM,min(0,c-ratio*d2));
clims = [min(min(IM)),max(max(IM))];
figure
colormap(gray)
imagesc(X,Y,IM,clims)
axis equal
disp('If the display is not clear enough, change the value of *ratio*') 












