%% Darstellung von 3D Intensitaetsbildern
% 1.
load('CTHeadVolume');

% 2.
s = size(imageVolume);
montage(reshape(uint8(imageVolume), s(1), s(2), 1, s(3)));

% 3.
implay(reshape(uint8(imageVolume), s(1), s(2), 1, s(3)));


%% Segmentierung der Knochenstruktur
% 1.
imageVolumeBin = zeros(s(1), s(2), s(3), 'uint8');

% 2.
boneLevel = 200;

% 3.a)
idx = find(imageVolume > boneLevel);
% 3.b)
imageVolumeBin(ind2sub(size(imageVolumeBin), idx)) = 255;

% 4.
montage(reshape(imageVolumeBin, s(1), s(2), 1, s(3)));
% Qualitaet passt.


%% Nachbearbeitung der Segmentierung
%
imageProcBin = imageVolumeBin;
% 1.a)
sphere = strel('ball', 4, 4);
% 1.b)
imageProcBin = imdilate(imageProcBin, sphere);
% 1.c)
imageProcBin = imerode(imageProcBin, sphere);
figure
montage(reshape(imageProcBin, s(1), s(2), 1, s(3)));

% 2.a)
CC = bwconncomp(imageProcBin, 26);
% 2.b)
maxCompId = 1;
maxCompSize = numel( CC.PixelIdxList{maxCompId} );
for id = 2 : CC.NumObjects
    if numel( CC.PixelIdxList{id} ) > maxCompSize
        maxCompSize = numel( CC.PixelIdxList{id} );
        maxCompId = id;
    end
end
posBoneStructure = CC.PixelIdxList{maxCompId};
% 2.c)
imageProcBin = zeros(s(1), s(2), s(3), 'uint8');
imageProcBin(posBoneStructure) = 255;

% 3.
figure
montage(reshape(imageProcBin, s(1), s(2), 1, s(3)));

% 4.a)
imageProcBin = smooth3(imageProcBin);
% 4.b)
idx = find(imageProcBin > 120);
imageProcBin(ind2sub(s, idx)) = 255;

% 5.a)
CC = bwconncomp(imageProcBin, 26);
% 5.b)
maxCompId = 1;
maxCompSize = numel( CC.PixelIdxList{maxCompId} );
for id = 2 : CC.NumObjects
    if numel( CC.PixelIdxList{id} ) > maxCompSize
        maxCompSize = numel( CC.PixelIdxList{id} );
        maxCompId = id;
    end
end
posBoneStructure = CC.PixelIdxList{maxCompId};
% 5.c)
imageProcBin = zeros(s(1), s(2), 1, s(3), 'uint8');
imageProcBin(posBoneStructure) = 255;


%% Surface-Rendering der Knochenstruktur
% 1.
SFimageProcBin = isosurface(imageProcBin, 100);
SFimageVolumeBin = isosurface(imageVolumeBin, 100);

% 3.a)
imageVolume = smooth3(imageVolume);

% 2.
figure;
patchedVolume = patch(SFimageVolumeBin, 'FaceColor', 'blue', 'EdgeColor', 'none');
% 3.b)
isonormals(imageVolume, patchedVolume);
% 3.c)
view(3);
lighting phong;
camlight;
camlight(80, -150);
daspect([1 1 0.56]);

% 2.
figure;
patchedProc = patch(SFimageProcBin, 'FaceColor', 'blue', 'EdgeColor', 'none');
% 3.b)
isonormals(imageVolume, patchedProc);
% 3.c)
view(3);
lighting phong;
camlight;
camlight(80, -150);
daspect([1 1 0.5]);

% 4. Die Stufen sind zwar immer noch vorhanden, aber es sind erkennbar
% weniger.