% Movie cartoon of sinograms and back projections in x-ray tomography
% 34 slices through a skull and brain distributed across a 5 x 7 matrix,
% whereby the central element is first reserved for the rotating projection
% recorded during the scan, with a yellow bar moving down the image. As
% each new region is arrived at by the yellow bar, the border of the
% corresponding sinogram lights up. Once this is finished, the sinograms
% are replaced with the developing back projections.

% Number of sinogram images corresponding to different heights = 35
% Number of projections = 181

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

% Create an array of sinograms 7 across and 5 down. The middle one (4,3) is
% reserved for the original rotating transmission image, then the resulting
% layers of tomograms
width = 0.12; % width of sinogram/backprojection image
height = 0.17; % height of sinogram/backprojection image
stepcol = 0.14; % horizontal steps between columns
steprow = 0.2; % vertical steps between rows
%sinoarray = zeros(7,5,4);

% Generic basis for making a frame using 'fill'
xTh = 0.05;
xFr = [0 1 1 0 0 xTh 1-xTh 1-xTh xTh xTh];
yFr = [0 0 1 1 0 xTh xTh 1-xTh 1-xTh xTh];

LL = 2*180/37;

% Read in original slice files, modified to remove outer arcs at bottom of
% images and renamed from the public domain wikipedia source
% https://commons.wikimedia.org/wiki/Scrollable_computed_tomography_images_of_a_normal_brain_(case_2)

for slice = 1:34
    str1 = 'skull/Computedtomographyofhumanbrain';
    str2 = num2str(sprintf('%02d',slice));
    str3 = '.png';
    imgStr = [str1 str2 str3];
    temp = imread(imgStr);
    myIm(:,:,slice) = rgb2gray(temp);
    Radon(:,:,slice) = radon(myIm(:,:,slice));
end

% Read in transmission projection image files
for proj = 1:75
    str1 = 'skull/3dCTscananimation-';
    % Start at file 20, close to skull being face on
    imgIndx = 1+mod(18+proj,75);
    str2 = num2str(sprintf('%02d',imgIndx));
    str3 = '.tiff';
    imgStr = [str1 str2 str3];
    temp = imread(imgStr,'PixelRegion',{[3,172],[31,200]});
    myIm2(:,:,proj) = rgb2gray(temp(:,:,1:3));
end


for theta = -1:180
    slice = 1;
    for row = 2:-1:-2
        for col = -3:3
            blx = 0.01+stepcol*(col+3)+(stepcol-width)/2; % bottom left x-coord
            bly = steprow*(row+2)+(steprow-height)/2; % bottom left y-coord
            pos = [blx bly width height];
            hold off
            subplot('position',pos);
            newplot
            R = radon(myIm(:,:,35-slice),0:theta);
            a = size(R,1);
            sinogramDummy = zeros(a,181);
            if (theta >= 0) % Not first frame
                sinogramDummy(:,1:theta+1) = R;
                finalIm = imrotate(sinogramDummy,-90);
            end
            if (theta < 0) % First frame
                if (row^2+col^2 == 0)
                    imshow(myIm2(:,:,1))
                    caxis([0 max(max(myIm2(:,:,1)))]);
                else
                    imshow(sinogramDummy);
                end
            elseif (row^2+col^2 ~= 0) % Not the central image
                imshow(finalIm);
                caxis([0 max(max(finalIm))]);
                frameFlagSt = 11 + (slice-1)*5;
                frameFlagEnd = 10 + slice*5;
                if (theta >= frameFlagSt) && (theta <= frameFlagEnd)
                    hold on
                    frame = fill(xFr*730,yFr*181,'y','LineStyle','none');
                    alpha(frame,0.7);
                end
                slice = slice+1;
            else % The central image - only 38 projections between 0 and
                % 180 degrees, so I fill in the frames in between by
                % merging neighbouring images with different weightings.
                
                % Relationship between theta and xx, the progression of the
                % central image from straight on (image 1) to 180 deg
                % around (image 38)
                xx = 1 + theta*37/180;
                yy = xx - floor(xx); % Goes from 0 to 1 in between images
                mergeIm = myIm2(:,:,floor(xx))*(1-yy) + myIm2(:,:,ceil(xx))*yy;
                imshow(mergeIm);
                caxis([0 max(max(mergeIm))]);
                hold on
                xline = [1 170];
                yline = 10+[theta theta]/2;
                sliceHeight = plot(xline, yline,'y','LineWidth',4);
                sliceHeight.Color(4) = 0.7;
            end
            axis square
            colormap(bone)
            hold off
        end
    end
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
    
end

% Fade out sinograms to black and central projection image to white
for fadeOut = 1:-0.05:0
    slice = 1;
    for row = 2:-1:-2
        for col = -3:3
            blx = 0.01+stepcol*(col+3)+(stepcol-width)/2; % bottom left x-coord
            bly = steprow*(row+2)+(steprow-height)/2; % bottom left y-coord
            pos = [blx bly width height];
            hold off
            subplot('position',pos);
            newplot
            if (row^2+col^2 ~= 0)
                Dummy = Radon(:,:,35-slice);
                finalIm = imrotate(Dummy,-90);
                imshow(finalIm*fadeOut);
                caxis([0 max(max(finalIm))]);
                hold on
                slice = slice+1;
            else
                imshow(myIm2(:,:,38));
                caxis([0 max(max(myIm2(:,:,38)))]);
                hold on
                fill([0 171 171 0 0],[0 0 171 171 0],'w','LineStyle',...
                    'none','FaceAlpha',1-fadeOut);
            end
            axis square
            colormap(bone)
        end
    end
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end

pauseMov(vid,15);

for theta = 0:180
    slice = 1;
    for row = 2:-1:-2
        for col = -3:3
            if (row^2+col^2 ~= 0) % Not the central image
                blx = 0.01+stepcol*(col+3)+(stepcol-width)/2; % bottom left x-coord
                bly = steprow*(row+2)+(steprow-height)/2; % bottom left y-coord
                pos = [blx bly width height];
                hold off
                subplot('position',pos);
                newplot
                
                R = radon(myIm(:,:,35-slice),0:theta);
                backp = iradon(R,0:theta);
                imshow(backp);
                caxis([0 max(max(backp))]);
                axis square
                colormap(bone)
                slice = slice+1;
            else % The central image
                blx = 0.01+stepcol*(col+3)+(stepcol-width)/2; % bottom left x-coord
                bly = steprow*(row+2)+(steprow-height)/2; % bottom left y-coord
                pos = [blx-0.019 bly-0.019 width+0.038 height+0.034];
                hold off
                subplot('position',pos);
                newplot
                set(gca, 'Projection','perspective');
                x = [-1 1;-1 1];
                y = [1 1;-1 -1];
                z = [0 0;0 0];
                slice2 = ceil((theta+1)*34/181);
                for ii = 1:slice2
                    dumIm1 = radon(myIm(:,:,ii),0:180);
                    dumIm2 = iradon(dumIm1,0:180);
                    dumIm3 = dumIm2;
                    dumIm3(dumIm2 <=50) = NaN;
                    hhh = surf(x,y,z+ii*10,dumIm3,'FaceColor','texturemap',...
                        'EdgeColor','none');
                    rotate(hhh,[0 0 1],theta*2+144.1,[0,0,0]);
                    hold on
                    hhh.AlphaData =(dumIm3>40);
                    hhh.FaceAlpha = 'texturemap';
                    % Add a dummy line to stop resizing of viewed volume
                    % with skull rotation angle
                    dumLine = plot3([-1. 1.],[0 0],[0 0],'Color','w');
                    dumLine.Color(4) = 0.0; % alpha = 0; fully transparent
                end
                xlim([-1 1]);
                ylim([-1 1]);
                zlim([-10 350]);
                axis square
                axis off
                colormap(gray)
            end
        end
    end
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end
% Output the movie as an mpg file
close(vid);