Skeletonizing a 3d voxel data
显示 更早的评论
I have voxel data in a csv file. It has 3 columns representing x,y,z coordinates and 9000 rows. After csv read and scatter3, I am able to visualize the 3d object in xy and z coordinates. the task at hand is to convert it into a binary image and skeletonize the 3D plot (3D object) and get the bracnhes and the end points of the skeleton. I am an amateur and so any suggestions would be helpful!
采纳的回答
Rik
2019-8-20
I generated some random points to show how your code should work. The attached function is an edited version of a FEX entry. This function doesn't require any toolbox. If you have the image processing toolbox, you can use bwskel instead.
L=false(20,15,30);
x=randi(size(L,1),1000,1);
y=randi(size(L,2),numel(x),1);
z=randi(size(L,3),numel(x),1);
ind=sub2ind(size(L),x,y,z);
L(ind)=true;
skel=Skeleton3D(L);
6 个评论
Karthick Jayaraman
2019-8-20
Rik
2019-8-20
You can either provide the actual size, or use the code below to derive it from the coordinates.
L=false([max(x) max(y) max(z)]);
Karthick Jayaraman
2019-8-20
Rik
2019-8-20
For the other steps are not trivial. Not too long ago I had to segment the first two branches from a binary tree. The code below formed the basis for that. It is a very rough version which is not really ready for sharing and may need modification. So your mileage may vary...
function [branch,endpoint]=segment_tree_branch(tree,seed,old_branch_end)
%follow a tree branch in the 27-nhood until a junction is found, return the
%traversed voxels as a logical
[x,y,z]=ind2sub([3 3 3],[14 1:13 15:27]);
shifts=[x'-2 y'-2 z'-2];
T=tree;loopcount=0;
while true
loopcount=loopcount+1;
T(seed(1),seed(2),seed(3))=false;
x=seed(1)+shifts(:,1);y=seed(2)+shifts(:,2);z=seed(3)+shifts(:,3);
remove_ind=...
(x<1 | y<1 |z<1) | ...
(x>size(T,1) | y>size(T,2) |z>size(T,3));
x(remove_ind)=[];y(remove_ind)=[];z(remove_ind)=[];
valid_shifts=shifts(~remove_ind,:);
inds=sub2ind(size(T),x,y,z);
if sum(T(inds))>1
%Multiple pixels are true, so this is the junction point. Unless
%this is the first loop, in which case the second one might be from
%the second branch. Try to follow the point that is the most
%co-linear. That means the area of the triangle formed should be as
%close to 0 as possible.
if loopcount==1
%formula: TriangleSurfaceArea=abs(cross(AB,AC))/2
S=zeros(sum(T(inds)),1);
AB=old_branch_end-seed;
n=find(T(inds)');
for k=1:numel(n)
AC=shifts(n(k),:);
S(k)=sqrt(sum(cross(AB,AC).^2))/2;
end
[min_area,k]=min(S); %#ok<ASGLU>
seed=seed+shifts(n(k),:);
else
break
end
elseif sum(T(inds))==0
%edge of tree reached
break
else
%select neighbor as seed for next iteration
seed=seed+valid_shifts(T(inds),:);
end
end
branch= tree & ~T;
endpoint=seed;
end
Karthick Jayaraman
2019-8-28
Rik
2019-8-28
The attached function is a slightly improved version of what I previously posted. Still not ready for the FEX, but you may find it usefull. This function returns the labeled tree. It should be possible to determine the children branches from this.
I have the feeling you're trying to use this for the same thing as I am: segmenting the bronchial tree.
I can't really help you with calculating the area and volume of segments, because that means you need to undo the skeletonization. For my goal it is enough to have the centerline, and I don't even need the generation number (apart from the trachea and the main bronchi), so for me this does the trick.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Image Processing Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
