--- a
+++ b/Thoracic Organs Segmentation code/RibCageApproximation.m
@@ -0,0 +1,166 @@
+function [FinalCurve]=RibCageApproximation(I,Thres);
+%Draw a  closed curve which approximates the cross
+%section of the rib cage in CT slices.
+%      [FinalCurve]=RibCageApproximation(I,Thres);
+%       Input parameters
+%       Img=> the CT slice
+%       Thres=> the desired threshold for the ribs (bones) extraction.
+
+Img=squeeze(I);
+
+%apply  threshold
+BinaryImage=Img>Thres;
+%perform morphological open operation
+BinaryImage = imopen(BinaryImage,strel([1 1;1 1]));
+
+%find the pixels of the rib cage.
+[x,y]=find(BinaryImage==1);
+
+%define the center of the image.
+ImgCenter =round( [size(BinaryImage)]./2);
+
+
+%Estimate the vector between bone pixels and image centre
+vector = [x-ImgCenter(1),y-ImgCenter(2)];
+%Vector's angle in degrees.
+Angle = mod(360+180/pi*atan2(vector(:,2),vector(:,1)),360);
+
+
+%In each directions (0, 10,...., 360 degrees), define the bone pixel which is closer to the image
+%center and keep it for the interpolation process.
+
+Theta=[0:10:360];k=0;
+
+Ind=cell([length(Theta)-1,1]);
+for i =1:length(Theta)-1
+    [Ind{i,1}] = find(Angle>=Theta(i) & Angle<Theta(i+1));
+    if ~isempty(Ind{i,1})
+      IndMinDist = FunctionMinDistance(x(Ind{i,1}),y(Ind{i,1}),ImgCenter(1),ImgCenter(2));
+      k=k+1;
+      IndOfMinDistance(k,1)=Ind{i,1}(IndMinDist(1),1);
+       
+    end 
+end
+
+
+%Set points with the minimum distance from the center of the ribs cavity in an order
+[xi, yi]=PixelsInOrder(x(IndOfMinDistance),y(IndOfMinDistance),ImgCenter(1),ImgCenter(2));
+
+%Cubic interpolation
+[xii, yii] = Curvefitting(xi,yi,0.2);
+xii = round(xii);yii=round(yii);
+
+%Draw the final curve
+FinalCurve=zeros(size(BinaryImage));
+
+for i =1:length(xii)
+    FinalCurve(xii(i),yii(i))=1;
+end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%
+function  IndOfMinDistance = FunctionMinDistance(I,J,c1,c2)
+%This function calculates the distance between each contour (I,J) point and
+%the centre of the contour
+%and finds the edge with the  minimum distance 
+%      IndOfMinDistance = FunctionMinDistance(I,J,c1,c2)
+
+%      Input Parameters
+%      [I,J] points coordinates
+%      [c1,c2] coordinates of the centre of the contour
+
+MinDistance = zeros(length(I),1);
+
+MinDistance(:) = sqrt((I(:)-c1).^2+(J(:)-c2).^2);
+IndOfMinDistance =find (MinDistance <= min (MinDistance));
+ 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+ function [Final_x, Final_y]=PixelsInOrder(x,y,c1,c2);
+
+%Define a cartesian system xy with center[c1,c2],
+%Set the points [x,y] in a circular order according to theirs angles 
+%       [Final_x, Final_y]=PixelsInOrder(x,y,c1,c2);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ c1=round(c1);c2=round(c2);
+
+ vector = [x-c1,y-c2];
+ %estimate the angle between the cartesian axis x and the points [x,y]
+ Theta = mod(360+180/pi*atan2(vector(:,2),vector(:,1)),360);
+ 
+ %sort the [x,y] points 
+ [ThetaSorted,Ind] = sort(Theta);
+ 
+ Final_x=zeros(length(x),1);
+ Final_y=zeros(length(y),1);
+ Final_x(1:length(x)) =x(Ind);
+ Final_y(1:length(y)) =y(Ind);
+ 
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ function [xi,yi] = Curvefitting(x,y,RES)
+% Interpolate(cubic) the points of the image in order to create a close curve
+%     [xi,yi] =Curvefitting(x,y,RES)
+%
+%      RES: resolution desired
+
+%     
+%      Modified Code by A. Koulouri
+%      
+%      Chenyang Xu and Jerry L. Prince, 3/8/96, 6/17/97
+%      Copyright (c) 1996-97 by Chenyang Xu and Jerry L. Prince
+%      Image Analysis and Communications Lab, Johns Hopkins University
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% convert to column vector
+x = x(:); y = y(:);
+
+N = length(x);  
+
+% make it a circular list since we are dealing with closed contour
+ x = [x;x(1)];
+ y = [y;y(1)];
+
+dx = x([2:N+1])- x(1:N);
+dy = y([2:N+1])- y(1:N);
+d = sqrt(dx.*dx+dy.*dy);  % compute the distance from previous node for point 2:N+1
+
+d = [0;d];   % point 1 to point 1 is 0 
+
+% now compute the arc length of all the points  from  the first point [x(1),
+% y(1)].
+% we use matrix multiply to achieve summing 
+M = length(d);
+d = (d'*uppertri(M,M))';
+
+% now ready to reparametrize the closed curve in terms of arc length
+maxd = d(M);
+
+if (maxd/RES<3)
+   error('RES too big compare to the length of original curve');
+end
+
+di = 0:RES:maxd;
+
+xi = interp1(d,x,di,'cubic');
+yi = interp1(d,y,di,'cubic');
+
+N = length(xi);
+
+if (maxd - di(length(di)) <RES/2)  % deal with end boundary condition
+   xi = xi(1:N-1);
+   yi = yi(1:N-1);
+end
+
+function X = uppertri(M,N)
+% UPPERTRI   Upper triagonal matrix 
+%            UPPER(M,N) is a M-by-N triagonal matrix
+
+[J,I] = meshgrid(1:M,1:N);
+X = (J>=I);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+