Cartoon of a tomography experimental setup

Matlab code 

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