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);