% read manual segmentation and generate a 3D plot

clc 

close all

clear all

addpath('functions');

%%

imset='training';

patient=1;

%% read direcotory of images

folder=['dcom/',imset,'/ED',num2str(patient),'/images'];

[IMAGES, xthickness, ythickness, zthickness] = gatherImages(folder);

 [x_max,y_max,z_max]=size(IMAGES);

% for k=1:z_max

%     subplot(3,5,k)

%     imagesc(IMAGES(:,:,k));

%     colormap(gray);

% end

% adjust the resolution so that it falls in range [0-max_res]

max_res=3;

norm_res = floor(max(IMAGES(:))/(max_res)) ; 

Im = (IMAGES./norm_res);

figure

for k1=1:size(Im,3)

I1=Im(:,:,k1);   

subplot(3,5,k1);

imshow(I1);

end  

%% read directory of contours

current_dir=pwd;

contour_dir=['dcom/',imset,'/ED',num2str(patient),'/contours/LV'];

cd(contour_dir);

dfile = dir('*.txt'); 

numFiles = size(dfile,1);

% load contours

for k=1:numFiles

    contours{k}=load (dfile(k,:).name);

end

cd(current_dir);

% display images and contours

for k=1:z_max

    figure (1)

    subplot(3,5,k)

    % show image

    I1=Im(:,:,k);

    imagesc(I1);

    colormap(gray);

    hold on

    % read contour

    C1=contours{k};   

    x1=floor(C1(:,1));

    y1=floor(C1(:,2));

    plot(x1,y1,'r');

    % creat segmentation mask

    figure (2)

    subplot(3,5,k)

    % note that row==y and column==x

    LV_seg(:,:,k)=accumarray([y1 x1],1,[y_max x_max]);    

    imagesc(LV_seg(:,:,k));

    colormap(gray)

    title('segmnetation mask')

end

%% find edges of segmentation contours

for k=1:size(Im,3)

    figure(3);

    subplot(3,5,k);

    %Find edges in intensity image

    edge_seg = edge(LV_seg(:,:,k)); 

    imshow(Im(:,:,k) -2*edge_seg); 

    title('edges')

end

%% show three dimensional patch

    figure

    phi = LV_seg(:,:,end:-1:1);

    p3 = patch(isosurface(phi, 0),  'FaceColor','red', 'EdgeColor','none'); 

    isonormals(phi,p3)

    zlabel('z','FontSize',10,'FontWeight','bold'); 

    box on ;  alpha(.6)

    xrange = []; 

    yrange = [];

    axis([1 size(phi,2) 1 size(phi,1) 1 size(phi,3)]) 

    set(gca,'YTick',yrange,'fontsize',8)

    set(gca,'XTick',xrange,'fontsize',8)

    set(gca,'ZTick',round(10:size(phi,3)/5:size(phi,3)), 'fontsize',8)

    hold on

    colormap(gray); 

    dimz=size(Im,3);

    for k1=1:3:dimz

    Img2D_CrossSection =  Im(:,:,k1); 

    xImage = [1 size(phi,2); 1 size(phi,2)]; %# The x data for the image corners

    yImage = [1 1 ; size(phi,1), size(phi,1)]; %# The y data for the image corners

    zImage = k1 * ones(2,2);   %# The z data for the image corners

    surf(xImage,yImage,zImage,...    %# Plot the surface

     'CData',Img2D_CrossSection ,...

     'FaceColor','texturemap'); 

    end    

    % viewpoint spec: view(az,el)

    view(-30,13)

axis tight

camlight

camlight(-80,-10)

lighting gouraud