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!

[darova]

Can you attach the data?

[Karthick Jayaraman]

Hi, I dont have my file that I need to work but here is an example. in Book1.csv I used the following command for a 3D graph

air = dlmread('Book1.csv',' ');

scatter3(air(:,1),air(:,2),air(:,3))

Answer 1

[Rik]

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);

[Karthick Jayaraman]

It maybe a dumb question but how is L determined if you already have x,y and z coordinates?! Like the one I shared a csv file before

[Rik]

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]

Ya, I get that now. However, doesnt this code generate a skeleton in 3D logical array but is there a way to identify the xyz coodrinates. My main goal of the skeletonizing process is to find the number of branches, where exactly the brnaches originate, what is the distance between each branch and finally the area of the entire surface from start of the branch tip till the end of the branch (in original 3D object before skeletonizing) Does that make any sense or am in confusing>?

[Rik]

You can use findND to convert the logical array back to a list of coordinates.

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]

The purpose of skeletonizing is to determine where the branch starts and find the area and volume of that entire segment corresponding to that branch in 3D. In that process, I notice there are many 3rd and 4th degree branches (like a branch that originates from a branch from the main stem).

Any suggestion on how to avoid the 3rd degree branch and calculate the volume of that not just the branch but the entire segment of that particular 3D structure would be helpful.

[Rik]

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.