Extract Skeleton From Image

extract_skeleton.m

distSkel=MIJ.getCurrentImage;


% transform to a list
[xx,yy,zz]=meshgrid(1:size(distSkel,1),1:size(distSkel,2),1:size(distSkel,3));

% find the non zero values inside objects
nonZeroPts=find(distSkel>0);

skelList=[xx(nonZeroPts) yy(nonZeroPts) zz(nonZeroPts)];

%% plot the points in 3D (if a 3d image)
plot3(skelList(:,1),skelList(:,2),skelList(:,3),'b.')


%% plot them in 2D
plot(skelList(:,1),skelList(:,2),'r.')

%% apply a threshold
skelPrune=distSkel>8;

% perform a component labeling to identify the connected and disconnected structures
skelCL = bwlabel(skelPrune);

% find the largest component
compVoxCount=hist(skelCL(skelCL>0),1:max(skelCL(:)));
[maxSize,argmax]=max(compVoxCount);

% show just these objects
mainSkel=skelCL==argmax;
imagesc(mainSkel);

%% find junctions
juncImage=nlfilter(double(mainSkel),[3 3],@(x) sum(x(:)));
ljuncImage=zeros(size(juncImage));
ljuncImage(juncImage==3)=2; % straight segments (3 voxels in neighborhood, in front, self and behind)
ljuncImage(juncImage<3)=1; % terminal segments
ljuncImage(juncImage>3)=3; % connections
ljuncImage(juncImage==0)=0; % set everything else to 0
imagesc(ljuncImage)