% Movie cartoon of a tomography experiment

clear all; close all
vid = VideoWriter('tomoExpt.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 15;
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;

[x,y,z] = sphere(270);
[a,b,c] = cylinder(1,100);
% Initialized blank sinogram image
sinogram = zeros(225,225,4,'uint8');

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

% Set up spherical bowl for lens surfaces
[X,Y,Z] = sphere(250) ;
X(Z>-0.8) = NaN ; Y(Z>-0.8) = NaN ; Z(Z>-0.8) = NaN;
rL = 1.2;
pos1 = [0.05 0.2 0.65 0.6];
pos2 = [0.71 0.2 0.24 0.6];

for ang = 0:2.4:357.6 % 150 frames, double the number of projections
    hold off
    subplot('position',pos1);
    newplot
    
    % Draw Tintin
    s = surf(-x,y,z*1.25,'EdgeColor','none','FaceColor',skin,...
        'FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.1); % Head
    rotate(s,[0 0 1],ang,[0 0 0]);
    hold on
    s = surf(0.5*a+0.1,b*0.5,c-1.5,'EdgeColor','none',...
        'FaceColor',skin,'SpecularStrength',0.1,'DiffuseStrength',1); % Neck
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.16*x+1.03),y*0.1,z*0.1,'EdgeColor','none',...
        'FaceColor',skin,'SpecularStrength',0.1,'DiffuseStrength',1); % Nose
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.25*x-0.2),y*0.1-0.98,z*0.3-0.1,'EdgeColor','none',...
        'FaceColor',skin,'SpecularStrength',0.1,'DiffuseStrength',1); % Right ear
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.25*x-0.2),y*0.1+0.98,z*0.3-0.1,'EdgeColor','none',...
        'FaceColor',skin,'SpecularStrength',0.1,'DiffuseStrength',1); % Left ear
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.05*x+0.94),y*0.07+0.3,z*0.11+0.16,'EdgeColor','none',...
        'FaceColor','k'); % Left eye
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.05*x+0.94),y*0.07-0.3,z*0.11+0.16,'EdgeColor','none',...
        'FaceColor','k'); % Right eye
    rotate(s,[0 0 1],ang,[0 0 0]);
    s = surf(-(0.08*x+0.86),y*0.16,z*0.17-0.5,'EdgeColor','none',...
        'FaceColor','r','DiffuseStrength',0.25); % Mouth
    rotate(s,[0 0 1],ang,[0 0 0]);
    % Quiff
    for jj = -3:3
        s = surf(-((0.6-abs(jj)/35)*x+0.35),y*0.11+0.04*jj,z*(0.28-abs(jj)/50)+0.91,...
            'EdgeColor','none','FaceColor',skin,'DiffuseStrength',1,...
            'SpecularStrength',0.1);
        rotate(s,[0 0 1],ang,[0 0 0]);
    end
    axis off
    axis equal
    xlim([-4.05 10]);
    ylim([-1.5 1.5]);
    zlim([-1.5 2.5]);
    light('Position',lp,'Style','infinite');
    
    % Incident x-rays denoted by arrow
    mArrow4([-4 0 0],[-1.4 0 0],'color','y','stemWidth',0.4,...
        'tipWidth',0.61,'facealpha',0.7);
    
    % Read file
    str1 = 'skull/3dCTscananimation-';
    % Start at file20, close to skull being face on
    imgIndx = 1+mod(19+floor((ang)/4.8),75);
    str2 = num2str(sprintf('%02d',imgIndx));
    str3 = '.tiff';
    imgStr = [str1 str2 str3];
    myIm = imread(imgStr,'PixelRegion',{[3,172],[31,200]});
    
    % Scintillator screen
    scTh = 0.1;
    scWi= 3;
    scHe= 3;
    scPos = 2.005;
    plotcube([scTh scWi scHe],[scPos -scWi/2 -scHe/2],0.5,'g',0);
    % including transmission image
    myIm2 = myIm(:,:,1:3);
    XX = [2 2; 2 2];
    YY = [1.5 -1.5; 1.5 -1.5];
    ZZ = [1.5 1.5;-1.5 -1.5];
    hhh = surf(XX,YY,ZZ,myIm2,'FaceColor','texturemap','EdgeColor','none');
    alpha(hhh,0.64);
    hold on
    
    view(-44,10);
    
    % Camera body
    % Main body
    plotcube([camSize camSize camSize],[s2cam -camSize/2 -camSize/2],1,[0.2 0.2 0.2],0);
    % Screw ring thread outer surface
    rCamAp = 0.9*camSize/2;
    [cx, cy, cz] = cylinder(rCamAp,100);
    camAp = surf(cx+s2cam,cy,cz*camSize/16,'FaceColor',[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(camAp,[0 1 0],-90,[s2cam,0,0]);
    % Screw ring thread inner surface
    [cx, cy, cz] = cylinder(0.9*rCamAp,100);
    camAp = surf(cx+s2cam,cy,cz*camSize/16,'FaceColor',[0.7 0.7 0.7],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(camAp,[0 1 0],-90,[s2cam,0,0]);
    % Front edge of thread
    th = 0:pi/180:2*pi;
    xa1 = rCamAp*sin(th);
    ya1 = rCamAp*cos(th);
    xa2 = 0.9*rCamAp*sin(th);
    ya2 = 0.9*rCamAp*cos(th);
    xa = [xa1 -xa2]+s2cam;
    ya = [ya1 ya2];
    za = 0.0*xa + camSize/16;
    ring = fill3(xa, ya, za,[0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(ring,[0 1 0],-90,[s2cam,0,0]);
    
    % Camera sensor
    % including image projected on sensor surface
    xSens = s2cam-sensThick;
    plotcube([sensThick camSize/2 camSize/2],...
        [s2cam-sensThick -camSize/4 -camSize/4],1,[0.5 0.5 0.5],0);
    XX = [1 1; 1 1]*(xSens)-0.005;
    YY = [1 -1; 1 -1]*camSize/4;
    ZZ = [1 1;-1 -1]*camSize/4;
    hhh = surf(XX,-YY,-ZZ,myIm2,'FaceColor','texturemap','EdgeColor','none');
    alpha(hhh,0.5);
    % Pulsing row of pixels highlighting where sinogram generated
    pulseAlp = 0.4+0.3*sin(2*pi*ang/90);
    strip = fill3([1 1 1 1 1]*xSens-0.01,...
        [-1 -1 1 1 -1]*camSize/4,...
        [-1 1 1 -1 -1]*camSize/125-camSize/20,...
        'y','LineStyle','none');
    alpha(strip,pulseAlp);
    
    % Lens
    % 1st convex surface
    lens1 = surf(2*X+4.4,2*Y,2*Z+1.6,'FaceColor',myBlue,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    alpha(lens1,0.5);
    rotate(lens1,[0 1 0],90,[4.4,0,0]);
    % 2nd convex surface
    lens2 = surf(2*X+4.5,2*Y,2*Z+1.6,'FaceColor',myBlue,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    alpha(lens2,0.5);
    rotate(lens2,[0 1 0],-90,[4.5,0,0]);
    % Outer rim of lens
    [lx, ly, lz] = cylinder(rL,100);
    camAp = surf(lx+4.4,ly,lz/10,'FaceColor',myBlue,...
        'FaceLighting','gouraud','DiffuseStrength',1,'LineStyle','none');
    alpha(camAp,0.5);
    rotate(camAp,[0 1 0],90,[4.4,0,0]);
    hold off
    
    % Sinogram through eyes and bridge of nose
    subplot('position',pos2);
    im2 = imresize(myIm2,75*3/170);
    ii = 1+floor(ang/4.8);
    sinogram(1+(ii-1)*3,:,1:3) = im2(sum(101:103)/3,:,1:3);
    sinogram(2+(ii-1)*3,:,1:3) = im2(sum(101:103)/3,:,1:3);
    sinogram(3+(ii-1)*3,:,1:3) = im2(sum(101:103)/3,:,1:3);
    imshow(sinogram(:,:,1:3));
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
    hold off
end
% Output the movie as an mpg file
close(vid);