% Movie cartoon of a scanning XRF experiment of the van Gogh painting 
% "Patch of Grass", revealing an earlier painting of a peasant woman's 
% head. Images taken from "Visualization of a Lost Painting by Vincent 
% van Gogh Using Synchrotron Radiation Based X-ray Fluorescence Elemental 
% Mapping", J. Dik et al., Anal. Chem. (2008) 80, 6436–6442.

clear; close all
vid = VideoWriter('patchOfGrass.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]); 

[a,b,c] = cylinder(1,100);

myBlue = [0.4 0.44 0.73];
lp = [-2 -1 1];
view(-44,10);
set(gca,'Projection','perspective');

% Read image files and reshape them as required 
myIm = imread('patchOfGrass.png'); % Main painting 
myIm2 = myIm(:,:,1:3);
myIm3 = imread('vermillion.png'); % Fluorescence map of mercury L line 
myIm4 = imread('naplesYellow.png'); % Fluorescence map of antimony K line 
myIm5 = imread('vermillionBW.png'); % B/W version 
myIm6 = imread('naplesYellowBW.png'); % B/W version 
% Generate reduced-sized images of fluorescent maps with dimensions equal
% to the steps taken in the x- and y-directions (row and col) 
myIm7 = double(imresize(myIm5,1/23))/256; % resized to 1104/23 = 48 x 1083/23 = 47 pixels
myIm8 = double(imresize(myIm6,1/23))/256; % resized to 1104/23 = 48 x 1083/23 = 47 pixels

% Plotting areas of three components of video 
pos1 = [0.0 -0.1 0.91 0.63]; % Experimental setup 
pos2 = [0.05 0.55 0.45 0.45]; % Naples Yellow fluorescence map 
pos3 = [0.5 0.55 0.45 0.45]; % Vermillion fluorescence map 

% Blank starting images for the step-by-step generation of the two
% fluorescent maps 
blankIm1 = zeros(size(myIm3,1),size(myIm3,2),size(myIm3,3),'uint8');
blankIm2 = zeros(size(myIm4,1),size(myIm4,2),size(myIm4,3),'uint8');

pogScale = 350; % Conversion factor of "patch of grass" in pixels to units of cartoon 
pogHe = size(myIm,1); % Height of original PoG painting 
pogWi = size(myIm,2); % Width of original PoG painting 
scantl = [245 522]; % top left pixel coords of scan region of "Patch of Grass"
scanbr = [839 1124]; % bottom right pixel coords of scan region of "Patch of Grass"
scanWidth = scanbr(1) - scantl(1); % Width of area of painting that was scanned 
scanHeight = scanbr(2) - scantl(2); % Height of area of painting that was scanned 
scanWidth = scanWidth/pogScale;
scanHeight = scanHeight/pogScale;

pogY0 = scantl(1)/pogScale; % Needed to position PoG in scanned region 
pogZ0 = -scantl(2)/pogScale; % Needed to position PoG in scanned region 
pogWi = pogWi/pogScale;
pogHe = pogHe/pogScale;

for col = 0:46 % 47 columns
    for row = 0:47 % 48 rows
        hold off
        subplot('position',pos1);
        newplot
        axis off
        axis equal
        axis tight
        light('Position',lp,'Style','infinite');
        
        % Original painting "patch of grass" 
        pogTh = 0.1; % Thickness of canvas 
        pogXpos = 7.001; % Position of canvas in x-direction (along incident-beam axis) 
        if (round(col/2) == col/2) % Move painting up
            pogZpos = pogZ0-row*scanHeight/47;
        else % Move painting down 
            pogZpos = pogZ0-scanHeight+row*scanHeight/47;
        end
        pogYpos = pogY0+col*scanWidth/47;
        
        % Canvas... 
        plotcube([pogTh -pogWi pogHe],[pogXpos pogYpos pogZpos],...
           1,[0.95 0.77 0.5],0);
        % ... including original painting "Patch of Grass" 
        hold on
        XX = [7 7; 7 7];
        YY = [0 -pogWi; 0 -pogWi]+pogYpos;
        ZZ = [pogHe pogHe; 0 0]+pogZpos;
        hhh = surf(XX,YY,ZZ,myIm2,'FaceColor','texturemap','EdgeColor','none');
        
        % -------------------------------------------------------------
        % Color plume according to how much Naples yellow and Vermillion
        % there is, using the ultra-low-res images of the fluorescent maps 
        % myIm7 and myIm8 
        if (round(col/2) == col/2) % moving up 
            yelFrac = round(100*myIm8(47-row+1,col+1));
            redFrac = round(100*myIm7(47-row+1,col+1));
        else % moving down 
            yelFrac = round(100*myIm8(row+1,col+1));
            redFrac = round(100*myIm7(row+1,col+1));
        end
        if (redFrac+yelFrac <= 100)
            blueFrac = (100 - (redFrac+yelFrac))/2; 
        else 
            blueFrac = 0; 
        end 
        yelStr = int2str(yelFrac); 
        redStr = int2str(redFrac); 
        blueStr = int2str(blueFrac); 
        colStr = [blueStr,'b',blueStr,'k',yelStr,'y',redStr,'r'];
        
        plumeCol = colormix(colStr); % Use https://ch.mathworks.com/matlabcentral/fileexchange/41595-colormix
        
        % -------------------------------------------------------------
        % Generate plume cloud 
        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.025+pogXpos,ycoordcone1/1.6-0.0375,...
            zcoordcone1/2,'FaceAlpha',0.32,'FaceColor',plumeCol,... 
            'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(xrfPlume,[3 2 0],-90,[pogXpos+0.025 -0.0375 0])
        % Dummy horizontal and vertical lines to stop image resizing
        plot3([7.101 7.101],[-2.35 3.05],[0 0],...
            'color','w', 'LineWidth', 0.1);
        plot3([7.101 7.101],[0 0],[-3.22 2.5],...
            '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 = pogXpos-FZPpos; % Horizontal distance FZP to painting 
        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.37,'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.37,'FaceColor','y','LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(focRad2,[0 1 0],90,[FZPpos,0,0]); 
        
        % -------------------------------------------------------------
        % Bright focal spot on sample 
        spot = 0:pi/45:2*pi;
        aa = 0.02*cos(spot);
        bb = 0.02*sin(spot);
        cc = 0.0*aa;
        fill3(cc+pogXpos-0.0025,bb,aa,'y','LineStyle','none','FaceAlpha',1,...
            'FaceLighting','gouraud','DiffuseStrength',1);
        fill3(cc*2+pogXpos-0.002,bb*2,aa*2,'y','LineStyle','none','FaceAlpha',0.5,...
            'FaceLighting','gouraud','DiffuseStrength',1);
        
        % -------------------------------------------------------------
        % Create detector
        detBody = surf(a/2+pogXpos,b/2,c+3.2,'FaceColor',[0.5 0.5 0.5],'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(detBody,[3 2 0],-90,[pogXpos,0,0]);
        detFaceTh = 0:pi/45:2*pi;
        aa = cos(detFaceTh);
        bb = sin(detFaceTh);
        cc = 0.0*aa;
        detFace = fill3(aa/2+pogXpos,bb/2,cc+3.2,myBlue,'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(detFace,[3 2 0],-90,[pogXpos,0,0]);
        view(-44,16);
        hold off
        
        % -------------------------------------------------------------
        % Plot developing fluorescence image of Naples Yellow
        subplot('position',pos2);
        if (round(col/2) == col/2) % Map developing vertically up 
            mapYpos = 47-row; % Starts at 47 and ends at 0
        else % Map developing vertically down 
            mapYpos = row; % Starts at 0 and ends at 47
        end
        mapXpos = col; % Starts at 0 and ends at 46 moving horizontally 
        % Define region of original fluorescence map that corresponds to a
        % single pixel in the ultra-low resolution images used for the
        % scanning 
        strow = 23*(mapYpos)+1;
        endrow = 23*(mapYpos+1);
        stcol = 23*(mapXpos)+1;
        endcol = 23*(mapXpos+1); 
        % Add this region to the initially black image showing the
        % progression of the fluorescence map 
        blankIm1(strow:endrow,stcol:endcol,1:3) = myIm4(strow:endrow,stcol:endcol,1:3);
        imshow(blankIm1); 
        
        % -------------------------------------------------------------
        % Plot developing fluorescence image of Vermillion
        subplot('position',pos3);
        if (round(col/2) == col/2) % Map developing vertically up 
            mapYpos = 47-row; % Starts at 47 and ends at 0
        else % Map developing vertically down 
            mapYpos = row; % Starts at 0 and ends at 47
        end
        mapXpos = col; % Starts at 0 and ends at 46 moving horizontally 
        % Define region of original fluorescence map that corresponds to a
        % single pixel in the ultra-low resolution images used for the
        % scanning 
        strow = 23*(mapYpos)+1;
        endrow = 23*(mapYpos+1);
        stcol = 23*(mapXpos)+1;
        endcol = 23*(mapXpos+1);
        % Add this region to the initially black image showing the
        % progression of the fluorescence map 
        blankIm2(strow:endrow,stcol:endcol,1:3) = myIm3(strow:endrow,stcol:endcol,1:3);
        imshow(blankIm2);

        frame = getframe(gcf);
        writeVideo(vid,frame);
        
%        hold off
    end
end
% Output the movie as an mpg file
close(vid);