% SIMART a Matlab program for the numerical SIMulation 
% of magnetic suspetibility ARTifacts in magnetic resonance
% imaging (MRI) - (c) S. Balac - 1999
% 
%  Stephane Balac 
%  Dept. de Mathematiques 
%  Univ. de Bretagne Occidentale 
%  29285 Brest - France
%  Stephane.Balac@univ-brest.fr
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESCRIPTION OF THE EXPERIMENT - these variables may be changed
%
% Description of the ball
% -----------------------
global r0; r0 = 1e-2;     % its radius in meter
global Xm; Xm = 3e-4;     % its magnetic susceptibility
%
% Description of the image sequence
% ---------------------------------
global  C;  C = [0 0 0];  % origin of the slice
global  c;  c = sqrt(sum(C.^2)); 
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;     % intensity of the gradients in tesla per meter
global ns; ns = [0 1 0];  % direction of the slice selection gradient
global nr; nr = [0 0 1];  % direction of the read-out gradient
global np; np = [1 0 0];  % 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 = 256;      % number of pixels (Np^2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning off % suppresses all subsequent warning messages
global beta; beta = (Xm*B0*r0^3)/3; 

NG=5;  % Quadrature formula for integration over the slice width
Gpts = [.906179, .538469, 0., -0.538469, -0.906179];
Gpds = [.236926, .478628, .568888, .478628, .236926];

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

times1=cputime;
IMG=[];       % array containing the image 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);
   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);
     IMG(ir,ip) = IMG(ir,ip) + Gpds(is)*Jac*BdB(3);
   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,IM)
axis equal

