skull sinograms and back projections matlab code

Back to imaging techniques animations

Download

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

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

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

frame = getframe(gcf);

writeVideo(vid,frame);

hold off

end

% Output the movie as an mpg file

close(vid);