This requires a File Exchange download,
load BW
dsk=@(r)strel('disk',r);
A=BW;
outer=~bwconvhull(A);
B=A|outer;
C=imopen(~B,dsk(2));
D=bwareafilt(C,1);
E=imfill(D|A,'holes'); %Closed shape
G=getInner(E,A); %Green (inner) boundary
R=bwmorph(imdilate(E,dsk(3)),'remove'); %Red (outer) boundary
rgb=255*cat(3,R,G,0*G) + (A&~(R|G));
imshow(rgb,[]); hold on
hold off
shg
function G=getInner(R,A)
[m,n]=size(A);
K=imdilate(A,ones(3))&imfill(R,'holes')&~A;
[I,J]=find(K);
[I,J]=orderPath(I,J);
V=[I,J];
for j=1:2
delta=vecnorm(diff(V,1,1),2,2);
v=find(delta>10);
v(end+1)=height(V);
discard=false(height(V),1);
for k=1:numel(v)
discard=trimfunc(V,v(k),discard);
end
V(discard,:)=[];
V=flipud(V);
end
G=poly2mask(V(:,2),V(:,1),m,n);
G=bwmorph( bwmorph(G,'thin',1) ,'remove');
end
function discard=trimfunc(V, i, discard)
T=V(i-40:i,:); %HARDCODE=40
t0=T(1,:);
d=T(end,:)-t0;
T=T-t0;
c=T/d;
dists=vecnorm(T-c.*d,2,2);
imax=find(dists==max(dists),1,'last');
discard(i-40+imax-3:i,:)=1;
end
function [I,J]=orderPath(I,J)
is=tspsearch([I,J],5);
I=I(is);
J=J(is);
k=find( hypot(diff(I), diff(J))>5 , 1);
I=circshift(I,-k);
J=circshift(J,-k);
end
