% Movie cartoon of a phase-contrast tomography experiment 

clear all; close all
vid = VideoWriter('phaseTomoExpt.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]);

% Read movie file 
vid2 = VideoReader('Edge_enhancement_with_distance_ZF_1080p.mov');
frames = read(vid2,[1 Inf]);

camSize = 2;

% Set up head and tail parts of zebrafish body as half ellipsoids 
[x,y,z] = sphere(270);
x(z<0) = NaN ; y(z<0) = NaN ; z(z<0) = NaN;

% Set up middle part of zebrafish body as modulated cylinder 
t = 0:pi/80:pi/2;
r = 0.2+0.8*(sin(t)).^2;
[a,b,c] = cylinder(r,100); 

% Set up eyes and otoliths 
[x2,y2,z2] = sphere(200);

% Set up backbone of zebrafish as modulated cylinder 
t2 = 0:pi/80:4*pi;
r2 = 0.2+0.8*(sin(2*t2)).^6;
[a2,b2,c2] = cylinder(r2,100); 

% Set up mirror faces 
xm1 = [0 0 0 0 0];
ym1 = [-1 1 1 -1 -1];
zm1 = [-1 -1 1 1 -1]; 


myBlue = [0.4 0.44 0.73];
skin = [0.945 0.761 0.49];
lp = [-2 -1 1];
view(-34,5);
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.02 0.15 0.75 0.7];
pos2 = [0.8 0.2 0.19 0.6];

% Apply the transform
for ZZ = 0:610 % 150 frames, double the number of projections  
    hold off 
    subplot('position',pos1);
    newplot
    set(gca, 'Projection','perspective');
    % Draw zebrafish embryo 
    surf(a/8,b/4,c-1,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4); % body part 1
    hold on
    surf(x/8,y/4,z*0.5,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4); % body part 2
    surf(x/40,y/20,-z*0.5-1,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4); % body part 3
    % Eyes
    surf(x2/10,y2/14+0.14,z2/7+0.25,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.25);
    surf(x2/10,y2/14-0.14,z2/7+0.25,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.25);
    surf(x2/25,y2/25+0.19,z2/25+0.25,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4);
    surf(x2/25,y2/25-0.19,z2/25+0.25,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4);
    % Otoliths
    s = surf(x2/35,y2/35-0.2,z2/25-0.05,'EdgeColor','none','FaceColor',[0.2 0.2 0.2],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.88);
    rotate(s,[1 0 0],29,[0,0,0]);
    s = surf(x2/35,y2/35+0.2,z2/25-0.05,'EdgeColor','none','FaceColor',[0.2 0.2 0.2],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.88);
    rotate(s,[1 0 0],-29,[0,0,0]);
    % Backbone 
    s = surf(a2/35,b2/35,c2-0.7,'EdgeColor','none','FaceColor',[0.5 0.5 0.5],...
        'FaceLighting','gouraud','DiffuseStrength',1,'FaceAlpha',0.4);

    axis off
    axis equal
    xlim([-4.05 7.1]);
    ylim([-1.5 1.5]);
    zlim([-1.5 4.1]);
    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);
    
    % Scintillator screen
    scTh = 0.1;
    scWi= 1.5;
    scHe= 1.5;
    scPos = 0.14+ZZ/305;
    plotcube([scTh scWi scHe],[scPos -scWi/2 -scHe/2],0.5,'g',0);
        
    view(-34,5);
        
    % Lens
    % 1st convex surface 
    lens1 = surf(2*X+2+ZZ/305,2*Y,2*Z+1.6,'FaceColor',myBlue,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);%,'FaceAlpha',0.35) ;
    alpha(lens1,0.5);
    rotate(lens1,[0 1 0],90,[2+ZZ/305,0,0]);
    % 2nd convex surface 
    lens2 = surf(2*X+2.1+ZZ/305,2*Y,2*Z+1.6,'FaceColor',myBlue,'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);%,'FaceAlpha',0.35) ;
    alpha(lens2,0.5);
    rotate(lens2,[0 1 0],-90,[2.1+ZZ/305,0,0]);
    % Outer rim of lens
    [lx, ly, lz] = cylinder(rL,100);
    camAp = surf(lx+2+ZZ/305,ly,lz/10,'FaceColor',myBlue,...
        'FaceLighting','gouraud','DiffuseStrength',1,'LineStyle','none');
    alpha(camAp,0.5);
    rotate(camAp,[0 1 0],90,[2+ZZ/305,0,0]);
    
    % 45-degree reflecting mirror 
    mface1 = fill3(xm1+4+ZZ/305,ym1,zm1,[0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',1);
    rotate(mface1,[0 1 0],45,[4+ZZ/305,0,0]);
    mface2 = fill3(xm1+4.1+ZZ/305,ym1,zm1-0.1,[0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',1);
    rotate(mface2,[0 1 0],45,[4.1+ZZ/305,0,-0.1]);
    mface3 = fill3(ym1*0.05*1.4142+0.05*1.4142+4+ZZ/305,xm1-1,zm1,...
        [0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',1);
    rotate(mface3,[0 1 0],45,[4+ZZ/305,-1,0]);
    mface4 = fill3(xm1+4+ZZ/305-1/sqrt(2),ym1,...
        zm1*0.05*sqrt(2)-1/sqrt(2)-0.1/sqrt(2),[0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',1);
    rotate(mface4,[0 1 0],-45,[4+ZZ/305-1/sqrt(2),0,-1/sqrt(2)]); 
    
    % Dot-dashed lines 
    plot3([-2 4+ZZ/305],[0 0],[0 0],'k','LineStyle','-.','LineWidth',1.6);
    plot3([4+ZZ/305 4+ZZ/305],[0 0],[0 2],'k','LineStyle','-.','LineWidth',1.6);

    % Camera body
    % Main body
    plotcube([camSize camSize camSize],[4+ZZ/305-camSize/2 -camSize/2 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+4+ZZ/305,cy,cz*camSize/16+2-camSize/16,'FaceColor',...
        [0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);    
    hold off 
    
    subplot('position',pos2);
    imshow(frames(:,:,1:3,ZZ+1))
        
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
    hold off
end
% Output the movie as an mpg file
close(vid);