% two-dimensional segmenter combined with Deep Learning
clear all   
close all
clc
addpath('functions')
%% load the image 
disp('Load MRI images');
imset='validation';
%imset='training';
%imset='online';
name1=['matFiles/',imset,'_dataES'];
load (name1)
%%
patient=15;
min_sn=sum(slice_per_patient(1:patient-1))+1;
max_sn=sum(slice_per_patient(1:patient-1))+slice_per_patient(patient);
slice_num=min_sn:max_sn;
%slice_num=max_sn;
%slice_num=round(mean(slice_num));
%slice_num=103;

% maximum iteration of active contour
max_its = 50;
intEweight=0.5;
DLweight=.25;
Dynamic_Window=0;

% enable plots
dis_ena=1;
save_ena=1;

% region of interest size
Mroi=70;

% get dicom info
dicom_path=['dcom/validation/ES',num2str(patient),'/images'];
%im_folder=dir(dicom_path1);
%dicom_path2=['dcom/validation/ED',num2str(patient),'/images'];
%full_dicom_filename = [dicom_path2 filesep im_folder(1).name];
para=get_dicominfo(dicom_path);


for k=1:length(slice_num)
display(['processing slice# ',num2str(slice_num(k))])
I=t_I(:,:,slice_num(k));
tic
if slice_num(k)>= max_sn
    Mroi=20;
end
%figure(1)
%maxplots=4;
%subplot(1,maxplots,1)
%imagesc(I);colormap(gray);title(['slice number=',num2str(slice_num)]);

%--- run DL to find ROI: LV location detection
% load parameters of Deep Learning for ROI
%load DL_ROI_params.mat;
%yROI=DLN(I,stackedAEOptTheta,inputSize,hiddenSizeL1,hiddenSizeL2,outputSize,netconfig);
%subplot(1,maxplots,2)
%imagesc(yROI);
%subplot(1,maxplots,2);hold on; plot(m_cnt(1),m_cnt(2),'r+')

% cut ROI from the image
%[subI,m_cnt]=mask2subImage(I,yROI,Mroi);
subI=t_Iroi(:,:,slice_num(k));

% changing the size of the ROI if we want, useful for apex slices
roi_x=round((size(subI,1)-Mroi)/2)+1:round((size(subI,1)+Mroi)/2);
roi_y=round((size(subI,2)-Mroi)/2)+1:round((size(subI,2)+Mroi)/2);
subI2=subI(roi_y,roi_x);
m_cnt=t_centers{slice_num(k)};
%subplot(1,maxplots,3)
%imagesc(subI);colormap(gray);hold on; plot(50,50,'r+')

% run DL-LV segmentation to get phi_0
file_name='DLconfigure/DL_LV_ES_params.mat';
load (file_name);
init_mask1=DLN(subI2,stackedAEOptTheta,inputSize,hiddenSizeL1,hiddenSizeL2,outputSize,netconfig);
init_mask2=clean_segs(init_mask1);
init_mask(:,:,slice_num(k))=remap_mask(init_mask2,size(subI)/2+1,subI);

% run segmentation
%figure(2)
%subplot(5,6,k)
[LV_seg1,phi,m_cnt] = region_seg_subPhi(I,subI2,m_cnt,init_mask2,max_its,intEweight,DLweight,Dynamic_Window,0,file_name);
LV_seg2=clean_segs(LV_seg1);
% resize mask to original size
out_mask = remap_mask(LV_seg2,size(subI)/2+1,subI);
LV_seg(:,:,slice_num(k))=out_mask;

toc
% compute metrics
manualPoints=t_contours{slice_num(k)};
[dm1(k),dm2(k),PD1(k),PD2(k),HD1(k),HD2(k)]=eval_metrics(I,m_cnt,LV_seg(:,:,slice_num(k)),manualPoints,para); 

% special care for slices with large Perp Distance
if PD2(k)>= 3 
    edited_seg=edit_prior_shape(subI2,LV_seg2,1);
    [dm11,dm22,PD11,PD22,HD11,HD22]=eval_metrics(I,m_cnt,edited_seg,manualPoints,para); 
    if PD22<PD2(k)
        dm1(k)=dm11;dm2(k)=dm22;PD1(k)=PD11;PD2(k)=PD22;HD1(k)=HD11;HD2(k)=HD22;
        LV_seg(:,:,slice_num(k))=remap_mask(edited_seg,size(subI)/2+1,subI); 
    end
%     [LV_seg1,phi,m_cnt] = region_seg_subPhi(I,subI,m_cnt,edited_seg,max_its,intEweight,DLweight,Dynamic_Window,0);    
%     LV_seg(:,:,slice_num(k))=clean_segs(LV_seg1);
%     [dm1(k),dm2(k),PD1(k),PD2(k),HD1(k),HD2(k)]=eval_metrics(I,m_cnt,LV_seg(:,:,slice_num(k)),manualPoints,para); 
end


end

%% evaluate segmentations 

% good contours
GP1=find(PD1<5);
GP2=find(PD2<5);

% averages 
AGP1=length(GP1)/length(PD1)
AGP2=length(GP2)/length(PD2)

PD1=PD1
PD2=PD2
average_PD1=mean(PD1(PD1<5))
average_PD2=mean(PD2(PD2<5))


DM1=dm1
DM2=dm2
average_dm1=mean(dm1(PD1<5))
average_dm2=mean(dm2(PD2<5))


% Hausdorff distance
average_HD1=mean(HD1(PD1<5))
average_HD2=mean(HD2(PD1<5))

[AGP1 average_dm1 average_PD1 average_HD1]
[AGP2 average_dm2 average_PD2 average_HD2]

%figure(5)
%subplot(1,2,1);stem(dm1);title('dice metric')
%subplot(1,2,2);stem(PD1);title('perpendicular distance')

%% display 2D segmentation
if dis_ena==1
    for k=1:length(slice_num)
    figure(1);subplot_tight(2,round(length(slice_num)/2),k)
    %imshow(subI2(:,:,slice_num(k)),'initialmagnification',200,'displayrange',[0 255]); hold on;
    imshow(t_Iroi(:,:,slice_num(k)),'initialmagnification',200,'displayrange',[0 255]); hold on;
    %contour(phi, [0 0], 'b','LineWidth',2);
    contour(LV_seg(:,:,slice_num(k)), [0 0], 'b','LineWidth',2);
    %contour(init_mask, [0 0], 'y','LineWidth',2);
    % contour(bwconvhull(init_mask), [0 0], 'r','LineWidth',2);
    %contour(t_yLV(roi_y,roi_x,slice_num(k)), [0 0], 'g--','LineWidth',2);
    contour(bwconvhull(LV_seg(:,:,slice_num(k))), [0 0], 'r','LineWidth',2);
    C1=t_contours{slice_num(k)};Ct=scaleContour(C1,t_centers{slice_num(k)},100);
    plot(Ct(:,1),Ct(:,2),'g--','LineWidth',2);
    %legend('auto seg','initial contour','manual seg','auto conv hul');
    init_mask_k=init_mask(:,:,slice_num(k));
    init_mask_r=remap_mask(init_mask_k,size(t_Iroi(:,:,slice_num(k)))/2+1,t_Iroi(:,:,slice_num(k)));
    hold on;contour(init_mask_r,[0 0],'y')

%     figure(2)
%     subplot_tight(2,ceil(length(slice_num)/2),k)
%     imagesc(t_I(:,:,slice_num(k)));colormap(gray);
%     hold on
%     contour(double(phiI(:,:,slice_num(k))), [0 0], 'r','LineWidth',2);
%     manualPoints=t_contours{slice_num(k)};
%     plot(manualPoints(:,1),manualPoints(:,2),'g','LineWidth',2)   
%     lvmaskI=remap_mask(t_yLV(:,:,slice_num(k)),t_centers{slice_num(k)},I);
%     contour(double(lvmaskI), [0 0], 'y','LineWidth',2);
%     contour(double(mask_rc(:,:,slice_num(k))), [0 0], 'b','LineWidth',2);
%     legend('auto','manual');
    end
 

%% 3d represenations
if length(slice_num)<2 || dis_ena==0
break
end

% V_LV_auto=LV_seg(:,:,(slice_num));
% V_LV_man=t_yLV(:,:,(slice_num));
% [x_max,y_max,z_max]=size(V_LV_man);

% find centers of masks
% for k=1:size(V_LV_man,3)
%     temp=regionprops(V_LV_man(:,:,k),'Centroid');
%     cnt_man(k,:)=temp.Centroid;
%     temp=regionprops(V_LV_auto(:,:,k),'Centroid');
%     cnt_auto(k,:)=temp.Centroid;
% end

% registeration to compensate mis-alignmnet 
% ref_man=cnt_man(round(z_max/2),:);
% for k=1:size(V_LV_man,3)
%     trans=fliplr(-cnt_man(k,:)+ref_man);
%     Vt_LV_man(:,:,k)=imtranslate(V_LV_man(:,:,k),trans);
% 
%     trans=fliplr(-cnt_auto(k,:)+ref_man);
%     Vt_LV_auto(:,:,k)=imtranslate(V_LV_auto(:,:,k),trans);
%     
% end


% create the bottom slice
% for k1=1:1
% Vt_LV_man(:,:,z_max+k1)=zeros(x_max,y_max);
% apex_cnt=round(ref_man);
% R=0;
% Vt_LV_man(apex_cnt(2)-R:apex_cnt(2)+R,apex_cnt(1)-R:apex_cnt(1)+R,z_max+k1)=1;
% end

% create the bottom slice
% for k1=1:1
% Vt_LV_auto(:,:,z_max+k1)=zeros(x_max,y_max);
% apex_cnt=round(ref_man);
% R=0;
% Vt_LV_auto(apex_cnt(2)-R:apex_cnt(2)+R,apex_cnt(1)-R:apex_cnt(1)+R,z_max+k1)=1;
% end

% top slice 
% top_cnt=round(ref_man);
% temp=zeros(x_max,y_max);
% temp(top_cnt(2),top_cnt(1),1)=1;
% temp(top_cnt(2)+1,top_cnt(1)+1,1)=1;
% Vt_LV_man2=Vt_LV_man;
% Vt_LV_man2(:,:,1)=temp;

% figure
% subplot(1,3,1)
% disp3d(flipdim(V_LV_man,3),0);
% title('original')
% axis off

% subplot(1,3,2)
% disp3d(flipdim(Vt_LV_man,3),0);
% title('translated')
% axis off

% subplot(1,3,3)
% disp3d(flipdim(Vt_LV_man2,3),0);
% title('translated')
% axis off

% subplot(1,3,3)
% disp3d(flipdim(Vt_LV_auto,3),0);
% title('translated')
% axis off 

%
% figure
% subplot(1,3,1)
% disp3d(flipdim(V_LV_man,3),9);
% title('original')
% axis off

% fz=7;
% subplot(1,3,2)
% disp3d(flipdim(Vt_LV_man,3),[fz,fz,fz]);
% title('translated')
% axis off
% view([-10 40]);

% fz=7;
% subplot(1,3,3)
% disp3d(flipdim(Vt_LV_man2,3),[fz,fz,fz]);
% title('translated')
% axis off
% view([-10 40]);

% fz=9;
% subplot(1,3,3)
% disp3d(flipdim(Vt_LV_auto,3),[fz,fz,fz]);
% title('translated')
% axis off
% view([-10 40]);

end
%%
%   fz=7;
%   figure  
%   disp3d(flipdim(Vt_LV_man,3),[fz,fz,fz]);
%   title('translated')
%   axis off
%   view([-10 40]);

% figure(1)
% subplot_tight(2,10,9)
% imshow(t_Iroi(:,:,slice_num(k)),'initialmagnification',200,'displayrange',[0 255]); hold on;
% contour(bwconvhull(LV_seg(:,:,slice_num(k))), [0 0], 'r','LineWidth',2);
% C1=t_contours{slice_num(k)};Ct=scaleContour(C1,t_centers{slice_num(k)},Mroi);
% plot(Ct(:,1),Ct(:,2),'g--','LineWidth',2);

%% show only one slice
%figure
%C1=t_contours{slice_num(k)};Ct=scaleContour(C1,t_centers{slice_num(k)},Mroi);
%showCurveAndPhi(t_Iroi(:,:,slice_num(k)), Ct,phi, slice_num(k))

%figure
%showCurveAndPhi(t_Iroi(:,:,slice_num(k)),Ct ,bwconvhull(LV_seg(:,:,slice_num(k))), slice_num(k))
%% save results
name1=['results/',imset,'/ES/case',num2str(patient)];
if save_ena==1
save (name1)
end