Creation of a set of sinograms and the subsequent tomogram slices via back projection

Matlab code 

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