%Before .m execution add paths
Path1= cd;  Path1=[Path1,'\OpticalFlow3D'];path(Path1,path);
Path1 = cd; Path1 = [Path1, '\MAP']; path(Path1, path);
Path1 = cd; Path1 = [Path1, '\TwoPhaseLevel Sets']; path(Path1, path);

%input images
%Inhale and Exhale  images
load InhaleImg.mat;load ExhaleImg.mat;

% InhaleImg=imresize(imread('InhaleImg.jpg'),0.3);
% ExhaleImg=imresize(imread('ExhaleImg.jpg'),0.3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Img=reshape([InhaleImg;ExhaleImg]',[size(ExhaleImg,2),size(InhaleImg,1),2]);

%Rib Cage Approximation
Thresh = Histogr(Img);


%MaxValue=max(max(max(Img)));
%Img = rescalingFunction(Img,0.16*MaxValue,MaxValue);
   
[RibCage]=RibCageApproximation(Img(:,:,1),Thresh);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Feature extraction 

%Optical Flow for the desplacement estimation
 i1= permute(squeeze(Img(:,:,1)),[1 2 3]);%Inhale Img
 i2= permute(squeeze(Img(:,:,2)),[1 2 3]);%Exhale Img
 Displacement = OpticalFlow3D(i1,i2); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum likelikehood pixel to belong Heart.  
[J,I] = MaxLikelihood(Img(:,:,1),RibCage,squeeze(Displacement));
 
 

InnerCurve = DrawCircle(4,I(1),J(1),size(Img,1),size(Img,2),1);

%Level Sets Process
[phi1,phi2]=TwoPhasesLevelSet(Img(:,:,1),RibCage,InnerCurve);

figure
 imshow(Img(:,:,1),[]); hold on
 c1 = contour(phi1,[0 0],'r','LineWidth',2.2);
 c2 = contour(phi2,[0 0],'g','LineWidth',2);
 title('Zero Level Set')
 hold off;

  
 [I_heart]=find((phi2)<=0);
 [I_lung]=find((phi1)<=0);

 SegmentOfThorax= zeros(size(Img));
 SegmentOfThorax(I_heart)=1;
 SegmentOfThorax(I_lung)=1;
 SegmentOfThorax= imfill(SegmentOfThorax);

 figure(1)
 imshow(Img(:,:,1),[]); hold on
 c = contour(SegmentOfThorax,[0 0],'m','LineWidth',2);
 hold off