Something like this might also work -
A=load('Image').Image;
T=findEdges(A,'top',10);
B=findEdges(A,'bottom',10);
C=findCenter(T,B);
imshow(A,[]);
hold on
for i=1:numel(T)
h=plot(T(i).x,T(i).y,'LineWidth',2);
plot(B(i).x,B(i).y,'LineWidth',2,'Color',h.Color);
plot(C(i).x,C(i).y,'--','LineWidth',2,'Color',h.Color);
end
hold off
function Line=findEdges(A,sense,n)
outside=imerode(~bwconvhull(A>150),strel('disk',5));
[~,S]=gradient(sgolayfilt(A,2,21,[],1));
if strcmp(sense,'bottom'), S=-S; end
B=(movmax(S,11,1)==S);
B(outside)=0;
B=imclose(B>0,strel('disk',4));
B=bwareafilt(B,n);
B=bwmorph(B,'thin');
rp=regionprops(B,'PixelList');
for i=numel(rp):-1:1
x=rp(i).PixelList(:,1);
y=rp(i).PixelList(:,2);
p=polyfit(x,y,4);
x=linspace(min(x),max(x),size(A,2));
y=polyval(p,x);
Line(i).x=x;
Line(i).y=y;
Line(i).ym=mean(y);
end
[~,reorder]=sort([Line.ym]);
Line=Line(reorder);
end
function Line=findCenter(T,B)
for i=numel(T):-1:1
Line(i).x=(T(i).x+B(i).x)/2;
Line(i).y=(T(i).y+B(i).y)/2;
Line(i).ym=[];
end
end
