%  Load all the DICOM data and make a video of the stack of data
%
%  Version 1.0, October 19, 2017, JAJ

%  Set overall parameters

dir_name='DICOM/ST00001/SE00001/';  %  Directory name
start_image=1;                      %  First image in series
end_image=747;                      %  Last image in series
frame_rate=50;                      %  Frame rate for playing back the movie

%  Set the dynamic range for the display

off_set=100;      %  Offset [Hu]
range=400;        %  Range to display [Hu]
map_values=128;   %  Number of gray level values

bone_off_set= -250;  %  Offset for showing the bones
bone_range=100;      %  Range fro showing the bones

%  Initialize figure

figure(1);
colormap(gray(map_values));
dicom_movie(end_image+1-start_image) = struct('cdata',[],'colormap',[]);

%  Read information for the first images

file_name='IM00001';
info=dicominfo([dir_name, file_name]);
dx=info.PixelSpacing(1);
dy=info.PixelSpacing(2);
Y = dicomread(info);
[Nx,Ny]=size(Y);

%  Make space for all the images

Y=zeros(Nx,Ny,end_image+1-start_image);
z_positions=zeros(end_image+1-start_image,1);

%  Loop through the images and read and display them

for i=start_image:end_image
    file_name=['IM00',num2str(floor(i/100)),num2str(floor(rem(i,100)/10)),num2str(rem(i,10))];
    info=dicominfo([dir_name, file_name]);
    Y(:,:,i) = dicomread(info);
    image(((1:Nx)-Nx/2)*dx,((1:Ny)-Ny/2)*dy,(double(Y(:,:,i))+off_set)/range*map_values)
    xlabel('Distance [mm]')
    ylabel('Distance [mm]')
    pos=sprintf('%5.1f, %5.1f, %5.1f', info.ImagePositionPatient(1), info.ImagePositionPatient(2), info.ImagePositionPatient(3));
    z_positions(i)= info.ImagePositionPatient(3);
    text(-150, -150, ['Image ', num2str(i),', Pos. (x,y,z) = ', pos,' mm'], 'Color', [1 1 1])
    axis('image')
    drawnow
    dicom_movie(i)=getframe;
end

%  Display the movie

movie(dicom_movie, 5, frame_rate);

%  Make a 3D surface from the data, flip the data upside down for right
%  view

disp('Making a volume of the data set')
data = smooth3(uint16(double(Y(:,:,end:-1:1))+bone_off_set)/range*map_values);

%  Display the data as an isosurface with end caps.

figure(3), clf
disp('Showing the data')
patch(isocaps(data,.5), 'FaceColor','interp','EdgeColor','none');
p1 = patch(isosurface(data,.5), 'FaceColor','blue','EdgeColor','none');
isonormals(data,p1);
view(3); 
axis vis3d tight
camlight left
colormap(gray(map_values));
lighting gouraud


