A full shape analysis for 3D images segmented in matlab in the format used by bwlabeln or watershed

ellipsoid_analysis.m

% this function (will) performs the ellipsoid based shape analysis on a 3D image
% the input is a labeled image probably from the bwlabel function
% the output is an array formatted like such
% labelId, volume, centerX, centerY, centerZ, extentsX, extentsY, extentsZ, pca1X, pca1Y, pca1Z, score1, score2, score3
function [out_table]=ellipsoid_analysis(in_image,out_file)
% create an array the same size as the image of x,y, and z points. This way we can calculate
% center of volume and other parameters by using multiplication
[xgrid,ygrid,zgrid]=meshgrid(1:size(in_image,1),1:size(in_image,2),1:size(in_image,3));

num_of_obj=max(in_image(:));

% assign the columns to the output table

col_names={'label','volume','cov_x','cov_y','cov_z','ext_x','ext_y','ext_z','pca1_x','pca1_y','pca1_z','score_1','score_2','score_3'};

num_of_outputs=length(col_names);
% initialize the output array
out_table=zeros(num_of_obj,num_of_outputs);

for cur_obj = 1 : num_of_obj
    % save the label and the volume
    out_table(cur_obj,getnameidx(col_names,'label'))=cur_obj;
    out_table(cur_obj,getnameidx(col_names,'volume'))=sum(in_image(:)==cur_obj);
    
    % get the x,y, and z positions of each voxel
    xpts=xgrid(in_image==cur_obj);
    ypts=ygrid(in_image==cur_obj);
    zpts=zgrid(in_image==cur_obj);
    
    out_table(cur_obj,getnameidx(col_names,'cov_x'))=mean(xpts);
    out_table(cur_obj,getnameidx(col_names,'cov_y'))=mean(ypts);
    out_table(cur_obj,getnameidx(col_names,'cov_z'))=mean(zpts);
    
    % extents analysis
    out_table(cur_obj,getnameidx(col_names,'ext_x'))=range(xpts);
    out_table(cur_obj,getnameidx(col_names,'ext_y'))=range(ypts);
    out_table(cur_obj,getnameidx(col_names,'ext_z'))=range(zpts);
    
    % shape tensor analysis
    covar_pts=cov([xpts(:) ypts(:) zpts(:)]);
    % in case it is only one point
    if (prod(size(covar_pts))<9)
        covar_pts=ones(3);
    end
 
    
    % save the first (largest) eigenvector as pca1, 
    % and the scores (eigenvalues) as score1...
    [evec,evals]=eig(covar_pts);
    evals=diag(evals); % scores are in a diagonal matrix
    first_comp=evec(:,3);
    out_table(cur_obj,getnameidx(col_names,'pca1_x'))=first_comp(1);
    out_table(cur_obj,getnameidx(col_names,'pca1_y'))=first_comp(2);
    out_table(cur_obj,getnameidx(col_names,'pca1_z'))=first_comp(3);
    
    out_table(cur_obj,getnameidx(col_names,'score_1'))=evals(3);
    out_table(cur_obj,getnameidx(col_names,'score_2'))=evals(2);
    out_table(cur_obj,getnameidx(col_names,'score_3'))=evals(1);
    
end
% make and write the header
header_txt=sprintf('%s,',col_names{:});
header_txt(end)=''; % delete the last comma
dlmwrite(out_file,header_txt,'');
% write the table
dlmwrite(out_file,out_table,'-append','delimiter',',');