Cartoon of XRF discovery of hidden van Gogh painting
Matlab code
% 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);