clc;clear all;close all;
%% Get all the images in the DICOM directory
dicom_path = 'CTdicom/';
files = dir([dicom_path '*.dcm']);
%% Read all the images
disp('start reading.');
clear DI;
for i = 1:length(files)
info = dicominfo([dicom_path files(i).name]);
% info has all the tags from the DICOM header
% InstanceNumber is the slice location
DI(:,:,info.InstanceNumber) = dicomread([dicom_path files(i).name]);
end
disp('end reading.') 
%%
%create mesh grids to make voxels
xlim = [size(DI,1), 1];
ylim = [size(DI,2), 1];
zlim = [size(DI,3), 1];
Resolution = 1;
ux = xlim(1) : -Resolution : xlim(2);
uy = ylim(1) : -Resolution : ylim(2);
uz = zlim(1) : -Resolution : zlim(2);
[X,Y,Z] = meshgrid( ux, uy, uz );
%%
%TODO hint: rescale DI values to [0,1]. Set a threshold for to extract the bone region.
maxval = 
minval = 
masks = 
%%
%TODO hint: use imerode/imdilate to clean the bone region
values = 
%%
%output thresholded result for visulization
mkdir('thresholded/');
for i=1:length(files)
imwrite(values(:,:,i),['thresholded/',num2str(i),'.jpg']);
end
%%
%TODO hint: look at all the thresholded results, try to find a bounding box to
%remove the background noise
cleanmask = 
values = values.*cleanmask;
%%
mkdir('cleaned/');
for i=1:length(files)
imwrite(values(:,:,i),['cleaned/',num2str(i),'.jpg']);
end
%%
figure;
ptch = patch( isosurface( X, Y, Z, values, 0.1 ) );
isonormals( X, Y, Z, values, ptch )
set( ptch, 'FaceColor', 'g', 'EdgeColor', 'none' );
set(gca,'DataAspectRatio',[1 1 1]);
xlabel('X');
ylabel('Y');
zlabel('Z');
view(180,0)
lighting( 'gouraud' )
camlight( 180,0 )
axis( 'tight' )
