coronary_refine.m~

function [] = coronary_refine( rpath )
% This function processes each probability map of coronary arteries under
% 'rpath'. The processing steps include but not limited to thresholding, 
% filling holes, thinning, detecting bifurcation or end points, 
% reconnecting disconnected branches, removing isolated branches, and 
% obtains a coronary artery tree finally.
% 
% Examples:
%   coronary_refine('path_of_parent_directory_containg_volumes')

% create output directory
wpath = fullfile(rpath, 'coronary');
if ~exist(wpath, 'dir'), mkdir(wpath); end

% processing each volume under rpath
img_list = dir(fullfile(rpath, '*.mha'));
for ii = 1:length(img_list)
    %% read mha volume
    img_path = fullfile(rpath, img_list(ii).name);
    [img_prop, img_info] = mha_read_volume(img_path);
    
    %% thresholdingCoro
    img_bin = img_prop >= (0.5*intmax('uint16'));
    % check the binary image obtained by considering it as a volume data,
    % and you can also store the binary volume into a single file (.mha file)
    w_info = img_info; % header information of volume to be written
    w_info.DataType = 'uchar'; % change the data type to uchar (uint8)
    mha_write(img_bin, w_info, 'path'); 
    
    %% filling holes
    % your code here ...
    img_bin = logical(img_bin);
    img_bin = bwmorph3(img_bin, 'fill');
    %% thinning
    sk = bwskel(img_bin); % skeleton
    % clear small pieces
    CC = bwconncomp(sk);
    numPixels = cellfun(@numel,CC.PixelIdxList);
    [~, idx] = sort(numPixels);
    for i = 1:length(idx)-30
        sk(CC.PixelIdxList{idx(i)}) = 0;
    end
    % plot the centerline of coronary artery by considering it as a set of
    % points, e.g. denote 'img_thin' as the result of thinning step
    sk0 = sk;
    
    %% detecting end points (bifurcation detection is in the next step)
    kernel = ones(3,3,3);
    kernel(2,2,2) = 100;
    edpt = convn(kernel, sk); % endpoint detection (conv)
    edpt = (edpt==101);
    edpt = edpt(2:size(edpt,1)-1, 2:size(edpt,2)-1, 2:size(edpt,3)-1);

    %% reconnecting disconnected branches & removing isolate points or branches
    sk = branchReconnect(sk, img_prop, 0.8, 45);
    CC = bwconncomp(sk);
    numPixels = cellfun(@numel,CC.PixelIdxList);
    [~, idx] = sort(numPixels);
    for i = 1:length(idx)-4
        sk(CC.PixelIdxList{idx(i)}) = 0;
    end
    vesselShow(sk0, sk);
    sk = branchReconnect(sk, img_prop, 0, 120);    
    v
    %% obtain coronary artery tree (by tracking or other methods)
    kernel(2,2,2) = 2; % convolution based bifurcation detection
    brch = convn(kernel, sk)>=5;
    brch = brch(2:size(edpt,1)-1, 2:size(edpt,2)-1, 2:size(edpt,3)-1);
    [r, c, v] = ind2sub(size(brch), find(brch));
    branches = [r, c, v];
    % search
    neighbor = [ 0,0,1; 0,1,0; 0,1,1; 1,0,0; 1,0,1; 1,1,0; 1,1,1; 0,0,-1;... 
        0,-1,0; 0,-1,1; 0,1,-1; 0,-1,-1; -1,0,0; -1,0,1; 1,0,-1; -1,0,-1;...
        -1,1,0; 1,-1,0; -1,-1,0; -1,1,1; 1,-1,1; 1,1,-1; 1,-1,-1; -1,1,-1;...
        -1,-1,1; -1,-1,-1];
    closed = brch;
    coro_tree = {};
    for j = 1:size(branches, 1)
        idx = branches(j, :);
        for i = 1:26 % pixels around joint
            cur = idx + neighbor(i, :);
            tempBranch = [];
            if sk(cur(1), cur(2), cur(3)) && ~closed(cur(1), cur(2), cur(3))
                closed(cur(1), cur(2), cur(3)) = true;
                tempBranch = [idx; cur];
                flag1 = 0;
                while 1 % branch found
                    flag = 0;
                    for k = 1:26
                        newnode = cur + neighbor(k, :);
                        if ~closed(newnode(1), newnode(2), newnode(3)) && sk(newnode(1), newnode(2), newnode(3)) &&...
                                ( isempty(find(brch(newnode(1)-1:newnode(1)+1,newnode(2)-1:newnode(2)+1,newnode(3)-1:newnode(3)+1),1)) )
                            flag1 = 1;
                            closed(newnode(1), newnode(2), newnode(3)) = true;
                            cur = newnode;
                            tempBranch = [tempBranch; cur];
                            flag = 1;
                        end
                    end
                    if ~flag
                        break
                    end
                end
            end
            if ~isempty(tempBranch)
                coro_tree = vertcat(coro_tree, tempBranch);
            end
        end
    end
    % 
    % plot the complete coronary artery tree in different color according 
    % to the ids of branches, e.g. denote 'coro_tree', a cell array, as the 
    % coronary artery tree obtained, and each element is a coordinate array
    % of a single branch
    coronary_show(coro_tree);
    clear;
%     %% save the tree obtained into a mat file (.mat)
%     coro_tree{1} = rand(10, 3); % for example, branch 1
%     coro_tree{2} = rand(12, 3); % for example, branch 2 ...
%     tree_name = split(img_list(ii).name, '.');
%     tree_name = [tree_name{1}, '.mat'];
%     tree_path = fullfile(wpath, tree_name);
%     save(tree_path, 'coro_tree');
end

end