The emergence of phase contrast with sample-detector distance

Matlab code 

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