Scanning x-ray fluorescence experiment cartoon
Matlab code
% Movie cartoon of a scanning XRF experiment
% Fossil and fluorescence-map images taken from "Phosphor imaging as a tool
% for in situ mapping of ppm levels of uranium and thorium in rocks and
% minerals", J.M. Cole et al., Chemical Geology vol 193, pp 127-136 (2003).
clear all; close all
vid = VideoWriter('xrfExpt.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
s2cam = 6.4;
camSize = 3;
sensThick = 0.16;
[xx,yy,zz] = sphere(200);
[a,b,c] = cylinder(1,100);
myBlue = [0.4 0.44 0.73];
lp = [-2 -1 1];
view(-44,10);
set(gca,'Projection','perspective');
xrf(1,:) = [12.97 1]; % Th L alpha1 energy [keV], intensity [arb. units]
xrf(2,:) = [16.2 0.8]; % Th L beta1
xrf(3,:) = [15.62 0.35]; % Th L beta2
xrf(4,:) = [3.69 1.4]; % Ca K alpha
xrf(5,:) = [4.01 0.47]; % Ca K beta
xrf(6,:) = [2.01 0.7]; % P K alpha
xrf(7,:) = [2.14 0.23]; % P K beta
xrf(8,:) = [13.61 0.4]; % U L alpha1
xrf(9,:) = [13.44 0.11]; % U L alpha 2
sigma = 0.14; % signal SD
phEnergy = 2:0.01:17; % Photon energy range of XRF signal
% Read files
myIm = imread('fishFossil.png'); % fossil fish
myIm2 = myIm(:,:,1:3);
myIm3 = imread('fishBWxrf.jpg'); % BW version of XRF yield
myIm4 = double(imresize(myIm3,0.02))/256; % scaled to 25 x 50 pixels
myIm6 = imread('fish.jpg'); % Fluorescence map
pos1 = [0.0 0.01 0.7 0.54];
pos2 = [0.05 0.6 0.4 0.34];
pos3 = [0.5 0.6 0.4 0.34];
blankIm = zeros(size(myIm6,1),size(myIm6,2),size(myIm6,3),'uint8');
for row = 1:25 % 50 sampling points wide
for col = 1:50 % 25 sampling points high
hold off
subplot('position',pos1);
newplot
axis off
axis equal
axis tight
light('Position',lp,'Style','infinite');
% Fossil fish
fossTh = 0.1;
fossWi= 2;
fossHe= 1;
fossPos = 7.001;
if (round(row/2) == row/2) % Move fossil to right
fossXpos = 1+1/50-(col)/25;
else % Move fossil to the left
fossXpos = -1-1/50+col/25;
end
fossYpos = (-12.5+row)/25-1/50;
plotcube([fossTh fossWi fossHe],[fossPos -fossWi/2+fossXpos -fossHe/2+fossYpos],...
1,[0.95 0.77 0.5],0);
% including fossil image
hold on
XX = [7 7; 7 7];
YY = [1 -1; 1 -1]+fossXpos;
ZZ = [0.5 0.5;-0.5 -0.5]+fossYpos;
hhh = surf(XX,YY,ZZ,myIm2,'FaceColor','texturemap','EdgeColor','none');
% Fluorescence plume
if (round(row/2) == row/2) % Map developing to the left
mapXpos = 50-col+1;
else % Map developing to the right
mapXpos = col;
end
myIm5 = imread('rainbow2.tiff','PixelRegion',{[1,512],...
[round(myIm4(row,mapXpos)*512/3),round(myIm4(row,mapXpos)*512)]});
conecoord1 = 0:pi/320:pi/2; % coordinates of light cone schematic
conecoord2 = 0:pi/150:2*pi; % coordinates of light cone schematic
[conecoord1,conecoord2]=meshgrid(conecoord1,conecoord2);
xcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
ycoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
zcoordcone1 = 7*(abs(cos(conecoord1))).^3;
xrfPlume = surf(xcoordcone1/1.6+0.04+fossPos,ycoordcone1/1.6-0.06,...
zcoordcone1/2,myIm5,'FaceAlpha',0.31+0.07*rand,'FaceColor',...
'texturemap','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(xrfPlume,[3 2 0],-90,[fossPos+0.04 -0.06 0])
% Dummy horizontal and vertical lines to stop image resizing
plot3([7.101 7.101],[0 -2],[0 0],...
'color','w', 'LineWidth', 0.1);
plot3([7.101 7.101],[0 0],[-1 1],...
'color','w', 'LineWidth', 0.1);
% -------------------------------------------------------------
% FZP
FZPpos = 5.5;
FZPth = 0:pi/100:2*pi;
imax = 14;
for jj = imax:-2:1
hold on
r1 = 0.125*(jj)^0.5;
r2 = 0.125*(jj-1)^0.5;
z_circle1 = r1*cos(FZPth);
y_circle1 = r1*sin(FZPth);
x_circle1 = 0.0*y_circle1;
z_circle2 = r2*cos(FZPth);
y_circle2 = r2*sin(FZPth);
x_circle2 = 0.0*y_circle2;
z_circle = [z_circle1 z_circle2];
y_circle = [y_circle1 y_circle2];
x_circle = [x_circle1 x_circle2];
zone = fill3(x_circle+FZPpos,y_circle,z_circle,'k','LineStyle','none');
end
% Create a cone of radiation focussing on sample
beamLength = fossPos-FZPpos; % Horizontal distance FZP to fossil
t = 0:0.05:beamLength;
r = (beamLength-t)/(beamLength);
[x,y,z] = cylinder(r,50);
focRad = surf(0.45*x+FZPpos,0.45*y,beamLength*z,...
'FaceAlpha',0.25,'FaceColor','y','LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(focRad,[0 1 0],90,[FZPpos,0,0]);
focRad2 = surf(0.45*x+FZPpos,0.45*y,-FZPpos*z,...
'FaceAlpha',0.25,'FaceColor','y','LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(focRad2,[0 1 0],90,[FZPpos,0,0]);
spot = 0:pi/45:2*pi;
aa = 0.02*cos(spot);
bb = 0.02*sin(spot);
cc = 0.0*aa;
fill3(cc+fossPos-0.0025,bb,aa,'y','LineStyle','none','FaceAlpha',1,...
'FaceLighting','gouraud','DiffuseStrength',1);
fill3(cc*2+fossPos-0.002,bb*2,aa*2,'y','LineStyle','none','FaceAlpha',0.5,...
'FaceLighting','gouraud','DiffuseStrength',1);
% Create detector
detBody = surf(a/2+fossPos,b/2,c+3.2,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(detBody,[3 2 0],-90,[fossPos,0,0]);
detFaceTh = 0:pi/45:2*pi;
aa = cos(detFaceTh);
bb = sin(detFaceTh);
cc = 0.0*aa;
detFace = fill3(aa/2+fossPos,bb/2,cc+3.2,myBlue,'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(detFace,[3 2 0],-90,[fossPos,0,0]);
view(-33,16);
hold off
% Fluorescence spectrum v scan
subplot('position',pos2);
flSpec = zeros(size(phEnergy)); % Intensity of fluorescence spectrum
for sig = 1:9
if (round(row/2) == row/2) % Map developing to the left
mapXpos = 50-col+1;
else % Map developing to the right
mapXpos = col;
end
if (sig>=1) && (sig<=3)
flSpec = flSpec + (xrf(sig,2)+myIm4(row,mapXpos))*...
exp(-(phEnergy-xrf(sig,1)).^2/(2*sigma^2));
else
flSpec = flSpec + xrf(sig,2)*exp(-(phEnergy-xrf(sig,1)).^2/(2*sigma^2));
end
sigout = poissrnd(1000*flSpec);
end
plot(phEnergy,sigout,'color','r', 'LineWidth', 2.0)
xlim([2 17]);
ylim([-150 2150]);
axl = get(gca,'XTickLabel');
set(gca,'XTickLabel',axl,'FontName','Helvetica','fontsize',18);
set(gca,'TickLength',[0.02, 1]);
set(gca,'linewidth',2);
xlabel('Photon energy [keV]');
ylabel('Fluorescence intensity [arb. units]');
hold off
% Plot developing fluorescence image
subplot('position',pos3);
if (round(row/2) == row/2) % Map developing to the left
mapXpos = 50-col+1;
else % Map developing to the right
mapXpos = col;
end
mapYpos = row;
stcol = 50*(mapXpos-1)+1;
endcol = 50*mapXpos;
strow = 50*(row-1)+1;
endrow = 50*row;
blankIm(strow:endrow,stcol:endcol,1:3) = myIm6(strow:endrow,stcol:endcol,1:3);
imshow(blankIm);
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end
end
% Output the movie as an mpg file
close(vid);