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